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We develop an analytic framework to understand fragmentation in turbulent, self-gravitating me- 
dia. In previous work, we showed how some properties of turbulence can be predicted by application of 
the excursion-set formalism. Here, we generalize this to understand fully time-dependent gravo-turbulent 
fragmentation and collapse. We show that turbulent systems are always gravitationally unstable in a prob- 
abilistic sense. The fragmentation mass spectrum, size/mass/density/linewidth relations of collapsing ob- 
jects, their correlation functions and clustering, the range of spatial scales over which fragmentation oc- 
curs, and the time-dependent rate of collapse/fragmentation (as a function of size/mass) are analytically 
predictable. We show how these depend on bulk properties of turbulence; fragmentation is promoted at 
higher Mach numbers and shallower power spectra. We also generalize the model to properly include 
rotation, complicated gas equations of state, collapsing/expanding backgrounds, magnetic fields, inter- 
mittency, and non-normal statistics (with inherently correlated fluctuations). This allows us to formally 
derive how fragmentation is suppressed with "stiffer" equations of state (e.g. higher polytropic index 7) or 
differently driven turbulence (solenoidal vs. compressive). The suppression appears at an "effective sonic 
scale" where bM(R s , p a n[Rs]) ~ 1, where p cl ; t is the (scale-dependent) critical density for fragmentation. 
Gas becomes stable against collapse below this scale for 7 > 4/3; however fragmentation still occurs 
on larger scales. We show that the scale-free nature of turbulence and gravity generically drives mass 
spectra and correlation functions towards universal shapes (observed in a wide variety of astrophysical 
phenomena), with weak residual dependence on many properties of the media. We find that correlated 
fluctuations on different scales, non-Gaussian density distributions, and intermittency have surprisingly 
small effects on the fragmentation process. We demonstrate that this is because fragmentation cascades 
on small scales are generically "frozen in" when large-scale fluctuations push the "parent" region above 
the collapse threshold; though they collapse, their statistics are only weakly modified by the collapse pro- 
cess. Finally, with thermal or turbulent support, structure develops "top-down" in time via a fragmentation 
cascade, but we show that significant rotational/angular momentum support reverses the sense of structure 
formation to "bottom-up" growth via mergers of bound clumps, and introduces a characteristic "maximal 
instability scale" distinct from the Toomre scale. 

Key words: hydrodynamics — instabilities — turbulence — star formation: general — galaxies: forma- 
tion — galaxies: evolution — protoplanetary discs — cosmology: theory 



1 INTRODUCTION 

Fragmentation and the collapse of gas under the influence of self- 
gravity in turbulent media is a process central to a wide range of as- 
trophysics. Galaxy formation within dark matter halos, the forma- 
tion of giant molecular clouds (GMCs) and structure within the in- 
terstellar medium (ISM), the formation of protostellar cores within 
GMCs, formation of binary and multiple stellar systems in proto- 
stellar disks, and planet formation within proto-planetary disks may 
all be fundamentally related to these basic physics (for reviews, 



see e.g. |Scal o & Elmegreen 2004; Elmegreen & Scalo 2004~| |Mac| 
|Low & Klessen|2004 1. As well, many processes as diverse as fu- 
sion within convective stars, grain growth in the ISM, and magnetic 
reconnection may be dramatically influenced by turbulent density 
and velocity fluctuations. 

Empirically, many independent lines of evidence suggest these 
phenomena are driven by some fundamental shared underlying 
physics and scale-free processes. For example, the mass functions 
(MFs) of stars at formation (IMF; for a review see Chabrier 2005 1, 
and of GMC s (|Blitz & Rosolowsky||2004 1, as well as protostel- 
lar cores (e.g. |Enoch et al.|2 008; Sadav oy"et al.|2010| >, star clusters 
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( |Portegies Zwart et al.|2010j>, and HI "holes" or underdense bubbles 
in the ISM ( |Oey & Clarke|1997||Walter & Brinks| 19991 and refer- 
ences therein) are nearly universal in shape. More remarkably, these 
mass functions are actually close to self-similar to one another, and 
even similar to the MF of (seemingly) completely unrelated sys- 
tems such as dark matter halos ! These all feature a power law-like 
range, with slopes close to dn/dM oc M~ 2 , i.e. the slope at which 
there is equal mass per logarithmic interval in mass (a generic ex- 
pectation of scale-free processes), with an exponential-like cutoff at 
low/high masses. Within these systems, quantities such as the mass- 
size relation (of GMCs, protostellar cores, ISM voids, massive star 
clusters, certain structures within turbulence regulating e.g. stel- 
lar temperature fluctuations, and dark matter halos) follow simple 
power-laws with near-universal slopes (see references above). The 
clustering (auto-correlation function) of young stars, proto-stellar 
cores, GMCs or dense (molecular) gas in the ISM, star clusters, and 
even galaxies also follows approximate power-laws (to first order), 
with - surprisingly again - apparently self-similar (near-universal) 
slopes (compare e.g. Zhang et al. 200 t]|Zehavi et all2 005 ; Stanke 
|et al.||2006| |Enoch et al.||2008| |Hennekemper et al.j|2008| |Kraus 
|& Hillenbrand|2008| |Scheepmaker et al.|2009| >. Extending this to 
the distribution of radial separations in binary and multiple stel- 
lar systems (where there is again a quasi-universal power-law dis- 
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tribution of separations) suggests that the same statistics may ex- 
tend self-similarly to fragmentation within proto-stellar disks (see 
Simon 1997, Kraus & Hillenbrand 2008, and references therein). 
And higher-order statistics (such as column density distribution 
shapes and velocity structure functions) of the turbulent ISM ap- 
pear consistent with simple scalings based on self-similar models 
that also apply to a wide range of laboratory MHD turbulence (e.g. 



Boldyrevetal.|2002||Ridge et al.|2006||Wong et al.|2008||Schmidt| 



et al.|2009[|Lombardi et al.|2010| >. Indeed, there is a long history 
of phenomenological models treating these phenomena as different 
aspects of fractal-like (i.e. strictly self-similar) systems (see e.g. 
Elmegreen 2002 , and references therein). 

Unfortunately, our theoretical understanding of the formation 
of self-gravitating structures within turbulent systems remains lim- 
ited. Early work on stability and fragmentation focused on smooth 



media dominated by thermal pressure and/or rotation (e.g. Jeans 



1902 



|Toomre|[T964l |Goldreich & Lynden-Bellj 1965; Li n et al 



1969 1. Subsequent analytic work extended this to include the ef- 
fect of turbulent "support" (energy, momentum, or ram pressure) 
in resisting collapse, in the dispersion relation for linear density 
perturbations in turbulent, rotating, and possibly magnetized disks 
( Chandrase khar|1951|[^ndervoort|1970l|Elmegreen|1987[[B^na7] 
zola et al. 1987 1. But these derivations still assumed that the media 
of interest were homogenous and steady-state, despite the fact that 
perhaps the most important inherent property of turbulent systems 
is their inhomogeneity. These analytic considerations provide no 
means to calculate the statistical properties of fluctuations in these 
turbulent media, let alone the statistics of complicated objects form- 
ing within those media or their time-dependence. 

This is not surprising: the systems of interest (fully de- 
veloped turbulence) are chaotic, non-linear, inhomogenous (with 
large stochastic fluctuations), span an enormous dynamic range 
(Reynolds numbers, or ratio of driving to dissipations scales as 
large as ~ 10 5 — 10 s ), intermittent (with shocks in super-sonic tur- 
bulence producing very large local perturbations), time-dependent, 
and include complicated thermal, magnetic, and radiative processes 
as well as differential rotation. As a result, most progress in the last 
couple decades has come from numerical simulations. 

Such simulations have led to a number of important break- 
throughs. Our understanding of the basic properties of astrophysi- 
cal turbulence is rapidly improving, and it appears to obey at least 
some surprisingly simple scalings in the velocity fields, albeit with 
significant intermittency (Ossenkopf & Mac Low 2002; Federrath 
|et al.|2010| [Block et al.|2010||Bournaud et al.|2010| >. In the ideal 
case of isothermal, non self-gravitating turbulence, at least, it ap- 
pears that density distributions can be described (approximately) 
as log-normal, with a dispersion that scales in a simple predictive 



manner with the large-scale compressive Mach num ber ([Vazquez- 



|Semadeni|19 94 ; Pado an et al.[l997l|Scalo et al.|19 98 , Ostrik er et al 
|1999[ >. Turbulence and the induced density fluctuations evolve (or 
decay, depending on the driving) on a crossing time (Pan & Scanna- 
|pieco|2 010 and references therein). For the problems of interest in 
star formation, the relation of turbulent driving, (super-sonic) Mach 
number, and density fluctuations in ideal "driven boxes" to the col- 
lapse rate of small self-gravitating regions is increasingly well- 
mapped ( Vazquez-Semadeni et al. 2003, Li et al. 2004; Padoan 
|& Nordlund||2011||Federrath & Klessen||2012| >. And in differen- 
tially rotating disks, some improvements to criteria for "effective 
stability" against fragmentation (beyond the original Toomre crite- 
rion) have been d eveloped ([G ammie 2001 ; Cai et al. 2008 ; H opkins | 
|& Quataert|2011 ; Elmegree n|2011) . Larger-scale simulations have 
attempted to follow the full "fragmentation cascade" in the ISM 



from GMCs to protostellar cores, with self-consistently driven tur- 
bulence from stellar feedback (e.g. Kim et al. 2002, Tasker 2011; 
Hopkins et al. 2012b a). These simulations and others like them 
have, individually, been able to reproduce many of the specific 
observations described above (e.g. mass functions, size-mass rela- 
tions, or correlation functions of some of the quantities of interest; 
see references above and Klessen & Burkert 2000; Jappsen et al. 
[2003][Krumholz et al.|201 l||Hansen et al.|2012| >. ~ 

But there are many caveats to the simulations. Numerical res- 
olution is limited, often well below the tremendous dynamic range 
in spatial scale, mass, and Reynolds numbers over which these 
processes operate. The space of interesting parameters describing 
the turbulence, fluid properties, and backgrounds is enormous and 
only a very small fraction can be surveyed. Sampling extremely 
rare events requires running simulations for durations which are in- 
feasible. And even when simulations are possible, the chaotic and 
non-linear nature of the systems means that interpreting the results, 
let alone extracting the important physics and understanding how 
to extrapolate it to systems not simulated - critical to understand 
the links between diverse astrophysical phenomena - is immensely 
challenging. 

Therefore considerable analytic theory has also been devel- 
oped, much of which has used the basic results of these simula- 
tions to develop predictions for a wide range of scales and tur- 
bulent properties. Simple analytic derivations underpin our un- 
derstanding of the approximately lognormal character of turbu- 
lent density PDFs (e.g. |Passot & Vazquez-Semadeni|1998[|Nord-| 
lund & Padoan 1999), and Passot & Vazquez-Semadeni ( 1998) ex- 
tended this further to predict the non-lognormal character in non- 
isothermal flows. A range of "cascade models" have been devel- 
oped that appear to successfully describe the scale-by-scale char- 
acter of hierarchical velocity fluctuations in fully non-linear turbu- 
lence (for a review, see She & Zhang 2009 1, and these have been 
increasingly applied to compressible density fluctuations ( Boldyrev 
|et al.|2002[[Kowal et al.|2007l|Liu & Fang|2008^ Analytic mod- 
els for star formation have been developed, which use the above 
scalings to calculate the mass or volume fraction exceeding some 
threshold density where self-gravity becomes important, and in 
turn use this to estimate the stellar initial mass function (Padoan 
|& Nordlundl|2002l |Hennebelle & Chabrier||2008"l |Veltchev et al. 
|201 1| >, mass distribution of clusters ( (Kle ssen & Burkert||200~n > 
and/or the integrated star formation rate (Krumhol z & McKee| 
|20031|Hennebelle & Chabrier|2011) . 

Building on these results, Hopkins (2012b) (hereafter Paper 
I) recently showed that the excursion-set formalism could be ap- 
plied to calculate the statistics of bound objects in the density field 
of the turbulent ISM. This is a general mathematical formulation 
for random-field statistics, well known from cosmological applica- 
tions as a means to calculate halo mass functions and clustering in 
the "extended Press-Schechter" approach from Bond et al.| ( [T99T| ). 
This is one of the most powerful theoretical tools in the study of 
large scale structure and galaxy formation, and the foundation for 
our analytic understanding of quantities such as halo mass func- 
tions, clustering, mergers, and accretion histories (for a review, see 
|Zentner|2007[ >. The application to the ISM therefore represents a 
potential major breakthrough, providing a means to calculate many 
quantities analytically that normally would require numerical sim- 
ulations. 

In Paper I, we focused on the specific question of GMCs in 
the ISM, and considered the case of isothermal gas with an ex- 
actly lognormal density distribution. We used this to construct cer- 
tain statistics of the "first-crossing distribution": the statistics of 
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bound objects defined on the largest scales on which they are self- 
gravitating. And we found that the predicted mass function and 
correlation functions/clustering properties agree remarkably well 
with observations of GMCs on galactic scales (formalizing many 
of the earlier calculations in e.g. Klessen & Burkert 200TJ. In |Hop-| 
|kins|(2012c| > (Paper II), we extended the formalism (with identical 
assumptions) to the "last crossing distribution" - specifically, the 
mass function of bound objects defined on the smallest scales on 
which they remain self-gravitating but do not have self-gravitating 
sub-regions (i.e. are not fragmenting). We argued that these should 
be associated with proto-stellar cores, and in Paper II showed that 
the resulting core MF agrees well with canonical Milky Way core 
MF and (by extrapolation) stellar IMFs. This connects earlier the- 
oretical work in Padoan & Nordlund (2002) and Hennebelle & 
Chabrier ( 2008 1 to a fully self-consistent galactic-scale framework. 
And in Hopkins ( 201 1 ) (Paper III) and Hopkins ( 20 1 2d I we showed 
how to calculate other properties of these "last-crossings" such as 
their correlation functions and dependence on the turbulent power 
spectrum. In addition to providing critical analytic insights into 
these quantities, this formulation allows us to simultaneously treat 
an enormous dynamic range in scales and consider extremely rare 
fluctuations, which are impossible to follow statistically in current 
simulations. However, the general applicability of these first mod- 
els is limited by some of the assumptions made; in Paper I-Paper III, 
we consider only isothermal, non-magnetized, isotropic gas, obey- 
ing strictly lognormal statistics with uncorrected fluctuations on 
different scales, and evaluate properties at a fixed instant in sta- 
tionary backgrounds (i.e. do not follow time-dependent collapse). 
We also considered only a narrow range of disk stability parame- 
ters, compressive-to-solenoidal ratios in the turbulence, and highly 
super-sonic Mach numbers. 

In this paper, we develop these initial models further, and gen- 
eralize them in many critical ways to enable the development of 
a theory of turbulent fragmentation applicable to a wide range of 
interesting astrophysical systems. Specifically, we generalize the 
models to include arbitrary turbulent power spectra, different de- 
grees of rotational support, complicated (multivariable) gas equa- 
tions of state, collapsing (or, in principle, expanding) backgrounds, 
magnetic fields and anisotropic media, intermittency, and non- 
Gaussian statistics (with or without inherently correlated fluctua- 
tions on different scales). We also develop a time-dependent ver- 
sion of the theory, both for statistically stationary and for glob- 
ally evolving backgrounds and collapsing self-gravitating "frag- 
ments." We show how, for all of these cases, quantities such as the 
mass spectrum of self-gravitating objects, their size-mass-density- 
linewidth relations, correlation functions/clustering, spatial/mass 
scales of the fragmentation cascade, and rate of formation, col- 
lapse, and fragmentation can be predicted. We present these as gen- 
eral models, applicable to a wide range of turbulent systems, but 
we also show that we can already reach a number of quite gen- 
eral conclusions regarding e.g. the "near-universality" of quantities 
like mass functions, size-mass relations, and correlation functions, 
as well as the conditions and characteristic scales (and range of 
scales) at which fragmentation occurs. 

1.1 Paper Overview 

The paper is organized as follows. A general illustration of the 
methodology is given in Fig. [T] and a number of important terms 
and variables that will be used throughout the paper are defined in 
Table[T| along with relevant equations for certain variables. 

In §[2}§[7] & §[9] we develop our methodology and gener- 
alize the theory in a number of important ways. § [2] reviews the 



basic methodology developed in Paper I-Paper III (for isothermal 
gas), shows how the model generalizes to rotating disks or collaps- 
ing subregions, and derives analytic solutions for mass functions 
and correlation functions. § |3|4| generalize this to polytropic and 
multi-variate equations of state. § [5] discusses how different driv- 
ing (compressive vs. solenoidal) enters the theory. § [^generalizes 
to include magnetic fields and anisotropy (in the velocity and/or 
density fields as well as collapsing structures). §[7] extends to in- 
termittent turbulence and highly non-Gaussian density statistics. In 
§ [9] we generalize all of these to be time-dependent, both follow- 
ing time-dependent collapse and evolution in statistically stationary 
background (§ |9.1[ > as well as time-dependent collapse and develop- 
ment of fragmentation within objects that are themselves collapsing 

(§E3- 

In § [8] we show and discuss the results, for fragmentation in 
fully-developed turbulence at fixed time. § |8.1| derives some basic 
key quantities and scaling relations, such as the "maximal insta- 
bility" and sonic scales, and mass-radius-density-velocity disper- 
sion relationships. § |8.2| considers the mass function of objects that 
form both on the largest self-gravitating scales and on the small- 
est scales in the fragmentation cascade, and how it depends on all 
of the parameters describing the medium (above). § |8.3| derives the 
distribution of the "dynamic range" of fragmentation, i.e. the extent 
of fragmentation cascades. § |8.4| presents the correlation functions 
and how they depend on properties of the medium. 

In § | 10| we consider the time-dependence of these results. 
§ |10.1|10.2| show how the global mass functions develop in time, 
how structures grow hierarchically and form fragments, and calcu- 
late the global rate of collapse of bound structures. 



1 1 0. 3 1 shows 

how sub-fragmentation develops in collapsing clouds with time, 
and how this depends on the initial properties of the cloud. § [TT| 
outlines the methodology needed to link these into a statistical en- 
semble of time-dependent "fragmentation trees." 

Finally, in §[T2]we summarize our results and discuss future 
work. Several more technical details are presented in the Appen- 
dices dAlGl. 



2 THE FIRST & LAST-CROSSING DISTRIBUTIONS: 
OUTLINE 

2.1 General Methodology 

In Paper I and Paper II we outline some of the methodology for ap- 
plying the excursion-set formalism to study the properties of self- 
gravitating structures in a turbulent medium (for simple cases). We 
briefly review this before generalizing this model to more compli- 
cated systems. 

If gas is isothermal, density fluctuations in both sub and super- 
sonic turbulence are (approximately) lognormal, so (by definition) 
the variable 5(x) = ln[p(x) / po] + 5/2, where p(x) is the density at 
a point x, po is the global mean density and 5 is the variance in In p, 
is normally distributed according to the PDFQ 

Po{8\S) = —L= expf-!^ (1) 



2S J 

More generally, we can evaluate the field 8(x\R), which is the S(x) 
field averaged around the point x with some window function of 
characteristic radius R. As shown in Paper I, this is also normally 
distributed, with a variance S(R) as a function of scale given by the 

1 The +5/2 term is just the subtracted mean in 5, required so that the in- 
tegral of pPg(p) correctly gives po with (8) = 0. If we consider the mass- 
weighted density distribution instead, then this becomes —S/2. 
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Figure 1. Illustration of the model and key concepts presented here. Left: The "barrier," or critical density above which a given region (with average density 
p(R) in a spherical radius R) inside a gaseous disk will collapse under self-gravity (Eqs. |13|14} . We label the dominant term resisting collapse at each scale: at 
large scales (> h) this is angular momentum, at intermediate scales h>R> R son j c this is turbulence, and at small scales R < i? son j c this is thermal and magnetic 
pressure. We label the sonic length S son ic (Eq. |45[ below which the Mach number M < 1, hence thermal support dominates, and the "maximal instability 
scale" (Eq. |44( near ~ h where fragmentation is most efficient. Center: A Monte-Carlo ensemble of several "trajectories." Each is an independent random 
realization of the density power spectrum, averaged in a window of variable size R about a random point in space x within the disk. Comparing to the barrier, 
most points in space (i.e. most of the volume) are not self-gravitating at any scale, but rare regions do cross the threshold. At large scales, density fluctuations 
must be suppressed by mass conservation, so p — > pq. At scales < R 80n j c , thermal pressure suppresses turbulent density fluctuations, so the smaller-scale 
fluctuations in each trajectory are damped and the curves become "frozen" at the values set by fluctuations over the scales where turbulence dominates. Right: 
Comparison of one high-density trajectory to the barrier. "First-crossing" is the largest scale R on which the region about the point x is self-gravitating, with 
enclosed mass given by Eq. |16| Multiple crossings below this scale represent independent sub-regions each collapsing within the parent, i.e. a fragmentation 
cascade. "Last-crossing" is the smallest scale on which the region remains self-gravitating - below this region no sub-scales are independently self-gravitating 
so there is no successive fragmentation. 



integral in Fourier space over the variance in fc-modes (i.e. the loga- 
rithmic density power spectrum) on scales < k < R~ l (discussed 
below) 

We define <S(x|/?) as the "trajectory" about a given random 
point in Eulerian space. If this is above some critical value B(R), 
it defines an "object of interest" on the relevant scale. In princi- 
ple, this can be almost anything; for this paper, we focus on self- 
gravitating collapse and fragmentation, so the obvious value of 
B(R) corresponds to the critical density above which a region on the 
scale R is self-gravitating. Of course, a given trajectory may cross 
B(R) many times in an arbitrarily narrow interval, corresponding to 
the region being self-gravitating on some scales, but not on others. 
We must therefore be careful what we mean by "self-gravitating" 
and how we select the relevant scales. 

In Paper I, we focus on the "first-crossing distribution": the 
distribution of self-gravitating regions defined on the largest scales 
on which they are self-gravitating. Physically, this corresponds to 
GMCs in the ISM (the largest self-gravitating gas structures). In 
Paper H, we define the "last-crossing" distribution instead, the dis- 
tribution of self-gravitating regions defined on the smallest self- 
gravitating scalesn 



2 Some subtleties related to the convolution of linear products of 
lognormally-distributed variables in real-space over a varying window (and 
important tests of the scale-by-scale assumption in numerical simulations) 
are discussed in Appendices |F|G| but these do not alter our results. 

3 Of course, once a region becomes self-gravitating and collapses, its den- 
sity can increase with time, departing from the lognormal density distribu- 
tion. We treat this in § 1 1 0-3 1 This does not, however, change our deriva- 



The properties of these regions can be derived in a Monte- 
Carlo manner. Recall, 5(x\R K .) is smoothed over some window - 
i.e. it is the convolution S(x \ R w ) = J dV W(\x' - xj, R w ) S(x'). So 
if we Fourier transform, we obtain <S(k|/?„.) = W(k|/?„.)5(k); the 
amplitude S(x | R w ) is simply the (window- weighted) integral of the 
contribution from all Fourier modes <5(k). Now begin at some suffi- 
ciently large scale Ri , where mass conservation implies Si (Ri) — > 
and r5i (x| Ri) — ?> 0. Then consider a "step" from this scale to R2 < 
Ri (S 2 > Si). We can then draw a new S 2 (R 2 1 <5i[i?i]) = S, + AS 
from the conditional PDF given AS = S 2 (R 2 ) — Si(Ri). This fol- 
lows from Fourier transforming the distribution of <5i, taking a 
"step" in Fourier space where we integrate the contribution from all 
contributing Fourier modes in the intermediate scales, and Fourier 
transforming back to real space (for a detailed discussion, see |Bond| 
|et al.|1991[|Zentner|2007| >. For the simplifying case where our win- 
dow function is a Fourier-space top-hatrlthis gives the simple result 

p(A+A*)dA*=-^e«p(-££)«A*) (2) 

Taking sufficiently small steps in AR to ensure convergence, the 
trajectory S(R) is then 

Rj>R t 

5(R,)=J2 AS '- ( 3 ) 



tions, since what we are interested in is the statistics of regions becoming 
self-gravitating via turbulent density fluctuations in the first place. 
4 Other window function choices are discussed in Appendix G 
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Table 1. Definitions of Important Terms & Variables Used Throughout This Paper 



Term Definition Eq. 

barrier critical density above which a region becomes self-gravitating (or any other criterion "of interest") 

trajectory random realization of the density field (density smoothed in successive scales R about a point x) 

sonic scale scale R below which the rms compressive turbulent velocities are sub-sonic (bv,(R) < c s ) 

maximal instability scale scale R at which fragmentation occurs most rapidly (near the disk scale height/turbulent driving scale) 

first-crossing largest scale on which a region is self-gravitating (trajectory crosses the barrier) 

last-crossing smallest scale on which a region is self-gravitating (has no self-gravitating/fragmenting sub-regions) 

crossing distribution the distribution of spatial scales (or corresponding masses) on whi ch first or last-crossings occur 

clustering strength (amplitude) of the two-point correlation function (Eq. 20| 



R spatial scale of a region (over which quantities are smoothed) - 
k Fourier mode k ~ \/R - 
p, p(R) density (p(R) is volume-averaged over a region of size R) 


S, S(R) variance in In (p) (volume-averaged within R) 
B, B(R) barrier (see above), for a region of size R 

v number of standard deviations corresponding to B (v = \B\/S 1 ^ 2 ) 


1 

4 


12 

n 


P probability distribution function (PDF) 
Pq Gaussian (normal) PDF 


1 




ft first-crossing distribution (defined above) 
fl last-crossing distribution (defined above) 


7 
4 




Pa\t(R) critical density for collapse (density associated with B(R)) 

po volume-averaged mean density of the system (mid-plane density in a disk) 


'i- 


h (exponential) disk scale height 

R c \ radius of an isolated cloud or turbulent box 


B 


re epicyclic frequency (re = 2flr g 1 d(r^O) / dr g where r g is the disk-centric radius, and re = re/Q) |l4| 
f! orbital frequency (Q = V c /r g where V c (r g ) is the circular velocity) 


Q Toomre Q parameter (defined on scale h, Q = (cr g [h] n)/(n GE gls )) 
Q 1 virial parameter for an isolated cloud/box (Q 1 = cr a [i? c i] 2 R d / (GM c j)) 


14 
14 




M, M(R) mass (M(R) is the critical mass for collapse within R, i.e. mass within R at p = p cr j t ) 

W(x, R) window function for smoothing the density field on scale R (W(k, R) is the Fourier transform) 


16 
12 




c s> c s (p, R) gas sound speed 

v,, v, (R) rms turbulent velocity dispersion (averaged on scale R) - 
i'a(p, R) Alfven speed (va =B/ •JA'K p) 

erg (p, R) total effective gas dispersion (erg = c\ + vj + v^) |l5| 


M Mach number M (R) = V, (R) /c s 

JV[ h Mach number at the disk scale height Mi, = v t [H]/c s 


b mean fraction of turbulent velocity in compressive (longitudinal) modes 
p turbulent spectral index over the inertial range (E(k) oc k~ p , vf <x R p ~ 1 ) 
7 gas polytropic index (c] oc p 7— ') 


l 2 

23 


J 


£(r | M) correlation function for objects of mass M, as a function of separation r 


2( 




ft intermittency parameter of the velocity field (/3 = 1 is non-intermittent, j3 = infinitely strong) 
13 p effective intermittency parameter for the density field (0 p = /3 1 / 3 ) 


35 




feoss(R) turbulent crossing time on a scale R ((cross (R) = R/vt(R)) 


^maximal, M maxima i maximal instability scale in radius or mass (M raaxima i = M(R maxilTm i)) 
Rsonio ^sonic sonic scale in radius or mass (M son j c = M(R son i c )) 


44 
45 





We can then construct an arbitrarily large Monte Carlo ensem- 
ble of trajectories, to evaluate various statistical properties of the 
mediumr] 



2.2 Analytic Solutions 



Zhang & Hui 1 2006 1 show that the first-crossing distributiorj^] fj 



for a normally distributed variable, as defined above, can be deter- 



Note that we have implicitly made an important assumption (in addition 
to assuming the PDFs are log-normal in the first place): namely, that the 
phases of different Fourier modes are un-correlated. This is unlikely to be 
true in detail in turbulent flows, where fluid parcels may have complicated 
correlation structures. In Appendix [F| we note that this cannot, in fact, be 
true if the PDF is lognormal on all scales; although we show that (based 
on simulations) it can be a good approximation for the quantities of interest 
if we focus on a portion (say, the high-density tail) of the PDF. We also 
show in § [3] that considering a non-isothermal gas forces us to explicitly 
include strong mode correlations (and discuss how to treat this). And in §|7] 
we consider some models for intermittency, which implicitly include an ad- 



justable degree of correlation structure, and therefore give some idea of how 
large the consequences may be. Even if the Fourier modes are intrinsically 
un-correlated, averaging in real-space itself introduces correlations between 
modes; this can be treated through the use of different window functions, 
discussed in AppendixjG] In any case, improving our understanding of the 
hierarchical structure in the density field is an important goal of ongoing 
and future numerical work, one which directly informs the class of analytic 
models we develop here. 

6 We define this, ffiS), as the fraction of trajectories which first cross 
above the critical density (i.e. become self-gravitating at S without being 
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mined purely analytically; and in Paper II we derive an analogous 
expression for the last-crossing distribution ft. This is given by the 
numerical solution to the Volterra integral equation: 



where 



f t (S) = gi(s)+ ds'Ms') g2 (s,s') 



Po(B(S)\S) 



gi(S)=\2 
8i(S, S') 



\dB\ 


5(5)1 


\dS\ 


5 J 



B(S)-B(S') + B(S) 



\dB\ 



5-5' 5 \dS\i 

P [B(S) - B(S') (5/5') | (5' - 5) (5/5')] 



(4) 

(5) 
(6) 



and 6(5) is the minimum value of the overdensity <5(x|6) which 
defines objects of interest (here, self-gravitating regions). The ex- 
pression for first-crossing is qualitatively similar but with some 
modifications (see Paper II): 

f f (S)=gi(S)+ f dS'f f (S')g 2 (S,S') (7) 
Jo 



where 



B{S)-B(S' 

5-5' ' " I dS 



+ 2 — h x 



(8) 
(9) 



6 [6(5)-5(5')|5'-5] 



These solutions are valid for any B and 5, provided that Po is 
Gaussian. We will use these expressions where possible, since they 
removes the need for Monte Carlo evaluation of trajectories, but 
emphasize that simply counting the fraction of trajectories which 
first cross B(S) (or fall below 6(5) for the last time) gives an iden- 
tical result (provided the number of trajectories is sufficiently large 
to converge). 

For the particularly simple case of a linear barrier (6 = Bo + 
fj,S), these yield the closed-form solutions: 



fe(S\B = B + nS) = 
f f (S\B = B + txS) = 




exp 



exp 



25 
If 
25 



(10) 
(11) 



Although the barrier is never exactly linear for any realistic turbu- 
lent regime, this is not a bad approximation in many regimes; for 
example, when d6/d5^>6/53>lwe recover the same limiting ex- 
pression for ft, and in the intermediate regime where d5/d5 <C 5/5 
we recover the same limiting expression for //. 

2.3 Dependence of the Variance and Critical Density on 
Turbulent Properties 

In Paper I we derive 5(6) and 6(5) from simple theoretical con- 
siderations for all scales in an isothermal, turbulent galactic disk. 
It is well-established that the contribution to density variance 
from the velocity variance goes as 5 ~ ln(l +M compressive), where 
-M compressive is the (compressive) Mach number (Federrath et al.| 



|2010| [Konstandi n et al.|2012b| and references therein). For a given 
turbulent power spectrum, this suggests that we can approximate 



so on any larger scales), per differential unit dS (i.e. per unit scale, with the 
scale variable being the variance). 



5(6) by summing the contribution from the velocity variance on all 
scales 6' > 6 



5(6): 



If 



(fc,6)| 2 ln[ 



1 + 



b 2 v?(k) 

d + K 2 k- 2 



dink (12) 



where W is the window function for the smoothing] v,(k) is the 
turbulent velocity dispersion averaged on a scale k (trivially related 
to the turbulent power spectrum), c s is the thermal sound speed, 
and b ~ 1 is the fraction of the turbulent velocity in compressive 
(longitudinal) motions. Here k is the epicyclic frequency; this must 
enter in the same way relative to c s as it appears in the dispersion 
relation for density perturbations (e.g. Toomre 19 77||Lau & Berlin] 
|1978| >; physically, this represents angular momentum suppressing 
large-scale fluctuations (and such a term must be present to ensure 
mass conservation, since the fluctuation amplitude must vanish suf- 
ficiently quickly on large scales). The constant b depends on certain 
properties of the turbulence, but we vary this below. For a deriva- 
tion of 5(6), see Paper I; we stress that the form we adopt is only an 
approximation, calibrated from numerical simulations, and discuss 
some alternatives in Appendix |E| 

Given S(R) = ln[p(6)/p ] + 5/2 defined above, then 6(6) 
follows from the dispersion relation for a density perturbation in 
a disk with self-gravity, turbulence, thermal pressure, and angular 
mom entum/shear (jyandervoort|1970| |Aoki et al.|1979| |Elmegreen] 
[19871 |Romeo|1992| >: 



5(6) = In 



Peril 



+ 



5(6) 



(13) 



where p C rit is the critical density above which a region is self- 
gravitating. This i^] 



Pcrit 
PO 



Q_ 

2k 



h\ yaJWh 2 R 
+ Rj laj(h) R + h 



(14) 



where po is the mean midplane density of the disk, h is the disk 
scale height, k = k/Q, — V2 (Q = V c /R) for a constant-circular 
velocity (V c ) disk, and Q = (a g [h] k) / (n GT, gll! ,) is the Toomre Q 
parameter. The dispersion <r g is 



■ c 2 s + {v}(R))+vl 



(15) 



(va is the Alfven speed). Again, the full derivation of these relations 
is presented in Paper I; for now, note that this ensures not only 
that a region is locally self-gravitating, but also that it can resist 
destruction by tidal shear (the k term) and/or energy input from 
shocks or the turbulent cascade (the v, term). 
The mapping between radius and mas^is 



M(R) \ 



t ,i\ Rl 



+ 



exp 



1 



(16) 



It is easy to see that on small scales, these scalings reduce to the 
Jeans criterion for a combination of thermal (c s ) and turbulent (v,) 
support, with M — (4tv/3) p C rit6 3 ; on large scales it becomes the 
Toomre criterion with M — 7rE C rit6 2 . 



7 For convenience we take this to be a i-space tophat inside k < 1 /R, which 
is implicit in our previous derivation, but we consider alternative window 
functions and their consequences in Appendix [g] 



Eqs. |14|16| are derived exactly for a disk with an exponential vertical 
profile. They are generically asymptotically exact at small and large \k\ and 
(comparing with numerical calculations) tend to be within ~ 10% of the 
exact solution at all \k\ for a wide range of observed vertical profiles (see 
|Kim et al.|2002) . 
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Recall ff(S) dS (or fe(S) dS) gives the differential fraction tra- 
jectories that have a first (last) crossing in a narrow range dS about 
the scale S[R] (corresponding to mass M = M[R]). Each trajectory 
randomly samples the Eulerian volume, so the differential num- 
ber of first (last)-crossing regions is related by V c \(M)dN(M) — 
V m fe(S[M])dS (where V c i(M) = M/p ait (M) is the cloud volume 
at the time of last-crossing and V to t is the total volume sampled). 
Hence, the mass function - the number density dn — dN /Vtot in a 
differential interval - is given by 



dn = Pcrit(M) I dS I 
dM M J{ 'IdMl 



(17) 



(for both ff and f£). 

The absolute units of the problem completely factor out; so the 
results generalize to broad classes of systems. However we need to 
specify the global stability parameter Q, as well as the basic proper- 
ties of the turbulence (e.g. the parameter b and the power spectrum). 
For convenience, we will generally focus on inertial-range turbu- 
lence E(k) oc k~ p , where p is the turbulent spectral index (usually 
p ~ 5/3 — 2)|^The normalization of this is defined by the Mach 
number at the scale h, M 2 , = (vf(h)) /c 2 s . 



2.4 Specification Within Collapsing "Regions" 

We will sometimes focus here on spherical, collapsing regions with 
radius R c i (defined generally, but analogous to GMCs or protostel- 
lar cores). In Paper III we derive the solution for the "two-barrier" 
problem for collapsing sub-regions within a parent distribution; for 
our purposes here, the important result is that if we consider the 
internal properties of initial clouds on some scale we can factor 
out the contribution from super-cloud scale fluctuations (effectively 
transforming to 5 — > S — So where So — S(R c i) without loss of gen- 
erality). We can therefore treat each region independently, with the 
variance "beginning" on the parent scale R c \ (i.e. treating this as a 
maximum radius as opposed to R — > oo); we do not need to know 
the behavior of the parent scales to model the continued evolution 
of the sub-region. 

If we consider small sub-regions, Ra <C h for inertial-range 
turbulence, we note that the variance and critical density simplify 
to 



S(R) 



R 

2 to Id 



din 



Pcrit ,l+M-a(R/Ra 
U l+M* 



po(Rc 



(18) 
(19) 



where Q' « 1 is trivially related to Q, but can for our purposes 
be taken to be an arbitrary virial parameter of the region, and 
po(Rd) = M c i/(47ri?|/3). The density ratio, then depends only on 
R /R c \ at fixed M c ] = M (R c \ ) . Thus, isothermal clouds form a one- 
parameter family in their behavior in M c i- 

2.5 Correlation Functions 

The auto-correlation function £ of a given population is defined as 
the excess probability of finding another member of the population 



9 As in Paper I, we note that this must turn over above R 2> h in a disk 
(to E(k) oak~ l ) to avoid an energy divergence, so we impose the flattening 
E(k) -> E(k) (1 + |i/2|- 2 )( l_ ' 7 )/ 2 , which gives a good fit to simulation 
results around the inertial scale (Bowman 1996 1. But as shown in Paper I, 
because the ft terms dominate Eq. |12|14| on these scales, the form of this 
turnover has weak effects on our results. 



within a differential volume at a radial separation r from one such 
member. 



1+£mm(7|M) ; 



(N{r\M)) 
(n(M))dV 



(20) 



where N(r\M) is the number with mass M in the separation r 
and n(M) is the average number density of objects of mass M. 
This is directly related to the conditional density PDF and first/last 
crossing distributions; specifically, to the probability that, given a 
first/last crossing on a given scale, the density field on a larger scale 
r has some value S, and then the probability that the field with some 
S on a larger scale contains multiple first/last crossings of the same 
mass. But these are both derivable from the existing field informa- 
tion - for example, by sampling the space of all possible "random 
walks" as described above. 

In Paper I-Paper III we derive the correlation function for both 
first and last-crossing distributions. For example, for last-crossings: 



l+&u(r|M) : 



fi(M\S ) 
,„ * MM) 

fe(M\S ^' 



MM) 



P*(<5o|S [r])d<5o 

Pb(S \S [r])dS (21) 



In the above, fe(M\So) is the last-crossing distribution (at mass 
M) determined for trajectories that begin with the initial condition 
So at scale S[r] = So- P*(So 1 5b M) is the probability of S(So) hav- 
ing the value So on the scale So, given that S(S[M]) = B(S[M]) 
- i.e. that there is a barrier crossing (a collapsing object) at the 
smaller scale, which is simply related to the conditional probabil- 
ity of S(S[M]) given the initial condition <5o(5o) by Bayes's theo- 
rem. Here fy (M \ So) is the conditional last-crossing distribution, i.e. 
the last-crossing distribution given an initial density So on a scale 
5oM; we show in Paper III that this is equivalent to fi(M) solved 
with the new initial condition for trajectories "beginning" at r with 
5^5-5 andB^B-<5 . 
For first-crossings, 



1+&m('-|^) = 



" W (ff(M\S ) 



//(AO 



q(So\Solr])dS (22) 



where the integral is over all So < 5 a i t (r), and q(So \ So) is a weight- 
ing factor defined in |Bond et al.| ( [T99T| ) as the probability that the 
overdensity at a random point, smoothed on a scale r, is So and does 
not exceed S C n t (r) on any larger smoothing scale. 

It is similarly straightforward to calculate the cross-correlation 
between objects of different masses[^]this amounts to replacing 
f(M\ So) 2 with /(Mi | So)f(M2 \ So). And the projected correlation 
function is just given by the line-of-sight integral (weighted by the 
appropriate number density). For details we refer to Paper I & Paper 
III. 

3 GENERALIZATION FOR POLYTROPIC GAS 

The above derivation considers only isothermal gas. We now gen- 
eralize this to locally polytropic gases: systems with c 2 oc p 7 ~ 1 . 

3.1 The Collapse Threshold 



Following Hennebelle & Chabrier (2009), it is straightforward to 
generalize our expression for the critical physical density above 

10 If fluctuations on different scales are uncorrelated; the correlation must 
be explicitly calculated from the Monte Carlo correlated random walk if 
there is explicit correlation structure in the Fourier modes of the density 
field. Wherever required, we use this method of calculation. 
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which a region will become self-gravitating, by replacing 

.7-1 



2 k 2/PV 



(23) 



h-^aJ (24) 



using the fact that c s (po) = co. This replacement in Eq,|15|makes 
Eq.[l4]an implicit equation for p cr j t : 

P^ = Q_ f, , M \ 4 (Pcrit/P0) 7 ~'+V, 2 (^) ft , _2* 

po ~ 2£^ + ^ I 4 + vj(h) 

which must in general be solved numerically. 

For collapsing spherical sub-regions (R c \ <C h) where we re- 
cover the Jeans criterion (i.e. drop the k terms), p cl -i t > k 2 [v,(k) 2 + 
Cj]/(4ttG), so Eq.fulsimplifies to 



Pcrit 
PO 



where we define Mid = Vt(R c i) /c s (po). It is straightforward to 
solve this numerically for p cr i t (if) . Note that so long as 7 < 2 there 
is always a unique solution. 

For Mc\{R/Rci) p ~ [ 1 turbulence dominates the sup- 
port and we obtain p cr j t « po(R/R c i) p ~ 3 , identical to the 
isothermal case (independent of 7). On the other hand 

for M 2 = Mlx{R/Rc\) p ~ l < 1 we obtain p crit » po(l + 
^2^-1/(2-7) ^ R ^- 2 /(z-y). the sub . sonic co ii a p Se density be- 
comes a steeper function of R as 7 increases. 
3.2 The Density PDF 

The generalization of the density PDF in the non-isothermal case is 
more delicate. Numerical studies have shown that it is no longer ex- 
actly lognormal in this case (see Passot & Vazquez-Semadeni 1998 
|Scalo et al.l[T998l [Nordlund & Padoan|19991|Ballesteros-Paredes| 
|et al.| |201 1 1. These authors point out, however, that the inviscid, 
unforced Navier-Stokes equations are invariant under the substitu- 
tion M 2 -> Ml s (p) = Ml (p/po)~ (7_1) (just the replacement of 
Cj(po) by Cj(p)). The same arguments driving the isothermal PDF 
to lognormality can be generalized, then, to predict the PDF of the 
form 



■[• 



P s = ln 



= exp 



25(7?, s] 



+ i>(R)s + <l>(R)) (26) 



where we replace the isothermal, constant variance S(R) = 
S(R, M[R]) with the local "effective variance" S(R, s) = 
S(R, A4zft[R, s]). The ip and (j> are simply determined for 
any S(R, s) by the normalization conditions J P(s) ds = 1 and 
j P(s) exp(s)ds = 1 (mass conservation). The replacement in S 
simply amounts to 



S(R, p) = 



\W(k, R)\ 2 ln \l + 



b 2 v 2 (k) 



C 2 (p[k]) + K 2 k- 2 



dln/fc (27) 



This should be valid so long as the turbulence obeys locality. In- 
deed, a variety of numerical experiments have shown that this is a 
good approximation, at least over the range 0.3 < 7 < 1.7 (see e.g. 
|Passot & Vazquez-Semadeni|1998||Scaio et al.|1998||Li et al.|2003> 

3.3 Constructing a Random Walk 

We now need to consider how we evaluate trajectories 5, if the PDF 
is non-Gaussian. There is a considerable literature on this topic for 



cosmological non-Gaussianity (see 


Matarrese et al.|2000; Afshordi 


|& Tolley|2008||Maggiore & Riotto 


2010b and references therein); 



however, most of the methods discussed therein are appropriate 
only for very small deviations from Gaussianity (the cosmological 
case of interest), whereas for 7 significantly different from unity 



we are confronted with entirely non-normal PDFs. However, we 
can take advantage of the fact that in a sufficiently small range of s, 
Eq.[26]is locally Gaussian-like. 

As in the isothermal case, when we consider a "step" from 
R 1 — > Rz = R 1 — AR and correspondingly S2(Ri) = si + As, we 
transform to Fourier space, integrate the contribution from all 
modes (according to our window function) and transform back. We 
simplify this considerably first by the use of a top-hat Fourier-space 
window function. In the limit of infinitesimal steps As, this reduces 
to the coupled stochastic differential equation 

As,- = (S0 + AR, s + As)) - (si(R, s)) 

+ s/S[R + AR, Jt + As,] - S[R, Si ] TZt (28) 

where IZi is an independent Gaussian random variable, and s; is the 
value for a single trajectory i. This is, to linear order 



As^Atf + 7^-Atf+-As 



As; 



A<*)| [+ 7^/AS| s 

P) 



(29) 

(30) 
(3D 



And the expansions are straightforward (albeit tedious) to higher 
orders. This can be solved numerically for each step of the tra- 
jectories in a Monte Carlo ensemble, populating P[s, R], if we re- 
strict our "step size" in AR to sufficiently small increments such 
that the expansions are valid to the chosen order (we caution that 
this requires some care). We solve the full equation with an itera- 
tive method to ensure the appropriate constraint is satisfied for all 
trajectories at all radii. Essentially, this amounts to locally near- 
Gaussian behavior, with an "effective" variance S e = S(R, s), and a 
shift/bias in the mean (s) because Eq.[26]is equivalent to a Gaussian 
with mean ip(R)S(R, s) if S(R, s) were constant; the d/ds terms in 
the denominator account to leading order for the locally-introduced 
skewness and kurtosis from the dependence of S(R , s) on s. 

We have explicitly tested the above (for both first and second- 
order expansions), with arbitrary c s (Ro | po) and v,(R) to generate 
S(R, p) in Eq.[27] to confirm that the Monte Carlo PDF that results 
does indeed agree well with the analytic Eq.|26]at each R (even far 
out in the tails of the distribution, where P < 10~ 7 ). We find that for 
the range of 7 of interest, a linear expansion is sufficiently accurate 
to recover all of the interesting behaviors provided the dependence 
of S(R, s) on s is sufficiently weak that | ^ + ^7= if 1*>\ < 1; 

otherwise a more exact solution is required. 

We stress that here (7 7^ 1) the steps are always correlated at 
some level, because the density s (itself determined by fluctuations 
on other scales) enters the equation to determine As. If the steps 
were not correlated, the central limit theorem would drive the PDF 
back to a Gaussian. An important question for study in numerical 
simulations in future work, is whether or not this correlation struc- 
ture approximated above is truly universal, and whether or not it 
implies different intermittency structures (since intermittency also 
implies correlated fluctuations) from the isothermal case. 

4 BAROTROPIC GASES & BIVARIATE GAS 
EQUATIONS OF STATE 

It is straightforward to see that our derivation for polytropic gases 
trivially generalizes to any barotrope c s = c,,(p). The polytropic 
index 7 has no unique role in the above other than defining the 
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barotrope, so is easily replaced by a more complex function (for ex- 
amples of this, see e.g. Jappscn et al. 2005, Henncbelle & Chabrier 



|2lX^|VeTtchev et al.|20lT 

Moreover, since we separate the various scales R in the above, 
it is also trivial to replace a barotrope with any bivariate function 
c s — > c s (p, R) of the density and radius. Or any quantity which is a 
function of the density and/or radius, such as the total mass inside 
a region (a product of p and R 3 ), the average turbulent velocity dis- 
persion within the region scale (using v,(R)), or the surface density 
of the region. 

In this paper, we simplify our analysis by restricting to poly- 
tropes. However, these bivariate distributions may be extremely rel- 
evant for a number of astrophysical cases: for example if the equa- 
tion of state is influenced by a uniform photo-ionizing background, 
or turbulent shocks. By simply replacing c s — > c s (p, R) in Eqs. 14 
& |27| the appropriate PDF, barrier, and method of following the 
Monte Carlo tree trivially follow. 

5 DEPENDENCE ON DRIVING MECHANISMS 



In this model (as in previous analytic and numerical work; Hen 



Inebelle & Chabrier 2 0091 |Padoan & Nordlund|[20TT1 [Federrath 



& Klessen 20121, the dependence of our results on the turbulent 
driving mechanisms, at otherwise fixed parameters, is entirely en- 
capsulated in two parameters. First, the turbulent power spectrum 
(freely varied below). Second, the parameter b, which represents 
the mean velocity component in compressive (longitudinal) as op- 
posed to solenoidal (transverse) motions. This enters the equations 
for S(R) since longitudinal motion is (by definition) what actually 
drives variation in p (see Konstandin et al. 2012b). If the turbu- 
lent driving is purely compressive, b = 1. For pure solenoidal driv- 
ing, b = 1/3 (a simple geometric consequence in isotropic, three- 
dimensional turbulence). For random driving with no "preferred" 
mode, b = 1/2 ( jPadoan & Nordlund|2002); for more detailed be- 
havior in intermediate cases, see Feder rath et al.|p010^ ; |Price et al.| 
pOTT) ; [Konstandin et al.1<2012b| l. 

6 MAGNETIC FIELDS & ANISOTROPIC COLLAPSE 

It is, in principle, also straightforward to generalize our results for 
magnetized media, with a few important caveats. A wide range of 
numerical experiments have shown that both sub and super-sonic 
turbulence in magnetized media also develop PDFs that obey many 
of the above scalings; for example, for 7=1 the PDF is also lognor- 
mal IQstriker et al.| 19991 |Klessen|20001|Kowal et al.|2007||Lemas^ 
|ter & Stone|2009[|Padoan & Nordlund|201 1^ " |Provided the shape 
of the PDF remains consistent with the class of barotropic solutions 
above, then there are a few additional ways magnetic fields may 
modify our conclusions. It is possible that MHD effects modify the 
turbulent velocity power spectra; however, for nearly all conditions 
of astrophysical interest which have been studied in numerical sim- 
ulations, the shape of the velocity power spectrum remains close to 
p ~ 5/3 (within the range p ~ 4/3 — 2 that we survey in this paper; 
see references above and Kritsuk et al.|20lT| >. 

The consequences of magnetic fields "resisting collapse" de- 
pend on the field geometry, which must be assumed. At one ex- 
treme, we can follow previous analytic approximations (e.g. |Kim| 
|et al.|2002| [Hennebell e & Chabrier|2009) and consider fields which 



' 1 In fact, as sho wn injH opkins ( 20 1 2a), from these simulations and those 
in |Molina et al.|(2012) , the log-normal approximation typically becomes 
more accurate with increasing magnetic field strength. 



are locally tangled on very small scales. In this limit, the field 
acts as an isotropic pressure term. In the criterion for collapse p cr j t 
(Eq.|14}, we simply take 



](p,R)^c 2 (p) + (v 2 (R)}+v 2 A (p,R) 



(32) 



where c 2 (p) allows for a non-isothermal equation-of-state and va is 



the Alfven speed v A = |B| /47rp = 



Likewise, in calculating 



the variance S(R) (Eqs.ll2|& 27 1 we add the magnetic pressure term 



resisting collapse, c 2 (p) — > c s [p] + v A (p, R). We allow the Afven 
speed to vary with density and radius - the variance with density is 
an effective "equation of state," while that with radius simply fol- 
lows the magnetic power spectrum. Having specified this, the mag- 
netic fields are mathematically degenerate with a bivariate e qua- 
tion of state (c 2 



^(p,R)[l+p- l (p,R)]) 



Kim et al. 



(2002 1 and 



|Kunz & Mouschovias (20101 compare this simple approximation 
with much more detailed numerical simulations and discuss where 
it may (and may not) apply; they suggest that where applicable, a 
near-constant va is also a reasonable approximation, further sim- 
plifying the relations above. 

At the opposite extreme, consider locally ordered magnetic 
fields B(x[U) = Bij B. Now, there is a preferred field axis, which 
has some important consequences. We can decompose the field and 
turbulence into components along each axis: 



2 . 2 



2 



(33) 



where v,._,- is the rms x-component of the turbulent field, and v fl x, 
v, || are the rms perpendicular/transverse (two-dimensional) and 
parallel/longitudinal turbulent components with respect to the B 
field. Compressions and expansions in this regime feel magnetic 
pressure perpendicular to the field lines, but not along the field 
lines. Therefore the differential term in the variance S(R) becomes 
approximately 



-(*) 



C}{p[k]) + K 2 k- 2 



+ ■ 



C*(p[k])+V 2 A (p,k) + Klk-l 
V?||(&) 



■ 2 ( P [k]) + K 2 k- 



(34) 



If we do not a priori know the magnetic field orientation, or it 
is changing in time (as is typical in turbulent fields), then we can 
further approximate the time-average by assuming isotropic turbu- 
lence, i.e. vf j_ R 
netized medium 



(2/3) vf and v, 11 ~ (1/3) v 2 . In a strongly, mag- 



Eq. 



34 



then simply becomes « (1/3) vf /(Cj + 
K z k~ z ); this is a pure geometric correction, as we assume the par- 
allel turbulent component is able to introduce compressions but the 
perpendicular components do not (similar to the scaling of b with 
solenoidal vs. compressive driving). While obviously approximate, 
this appears to capture the most important result (for our purposes) 
of MHD turbulence simulations, namely the dispersion-Mach num- 
ber relationF 2 ! 

12 Eq. [34] with the assumption of isotropic turbulence, in an ideal driven 
turbulent box, predicts a density-Mach number relation 

5f»ln[l+fo 2 A4 2 (l/3 + 2/3c 2 /(c 2 + v A ))] 

= In [1 + b 2 M 2 (2/3 + /3 A ) (2 + p A )] 

so the i'A-dependence can be absorbed into a "renormalized" (/3-dependent) 
b value. If we compare this to simulation s in|Lemaster & Stone|(2009} ; 
|Padoan & Nordlund|pOTT) ; |Molina et al.|j2012) , we find it gives a good 
approximation to the dispersion-Mach number relation and density power 
spectra (comparing as described in AppendixjE}, from field-strengths /3 A ~ 
20 — 10 — 3 surveyed therein (equivalent to an "renormalized" b — ¥ 0.58 b — 
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It is possible in principle (though outside the scope of this pa- 
per) to remove the isotropy assumption, and allow for the random 
walk in the density distribution to be inherently anisotropic. This 
amounts to considering the turbulent modes in Fourier space not 
as isotropic (i.e. functions of k alone) but as vectors k. The proce- 
dure for treating this is straightforward in the Monte Carlo method. 
Consider evaluating the random field that is the compressive com- 
ponent of the local Mach number (i.e. just the local velocity field). 
This decomposes into v t , x , etc., as above. So instead of evaluating 
one Gaussian random field variable as a function of scale for each 
trajectory, we would simply evaluate three (independent) variables 
(the three components M x , M y , A4 Z ). We associate with each a 
variance Sm, x (R), Sm.v(R), Sm, z (R); these are determined from the 
power spectrum just as in the isotropic case, but measured sepa- 
rately for each component (see an application of this approach for 
the angular momentum content of dark matter halos in |Sheth ~&\ 
|Tormen|[2002| l. For example, on large scales in a disk, it may be 
more accurate to similarly decompose vf into azimuthal, radial, and 
vertical modes. These can obey different power spectra (as seen in 
some simulations; Block et al. 2010, Bournaud et al. 2010), and 
similar to Eq.[34]have different support: the k term applies to radial 
modes, a similar term in f2 (and its derivatives) applies to azimuthal 
modes, and no angular momentum term resists vertical modes. In 
what follows, for simplicity and clarity, we always evaluate quanti- 
ties in spherical annuli, but in principle one could allow for triaxial 
collapse with this method. 

With that in mind, we see that the lowest-order corrections 
from magnetic fields - within the context of our model assumptions 
(which are admittedly highly simplified) - are mathematically iden- 
tical to changes in the freely-varied parameters. In other words, for 
magnetic field strength /3a in either limit described above, the for- 
mal approach and results are identical to those for a "pure hydrody- 
namical" case with appropriately re-mapped values of p, 7, a g , and 
b. For reasonable field strengths which have been explored numer- 
ically, the re-mapped "effective" values fall well within the range 
we treat as freely varied. Therefore, we will not explicitly denote 
magnetic vs. non-magnetic cases below, but note that the simple 
relations above can be used to determine the approximate effects 
for any specific case of interest (the general sense of increasing 
magnetic field strength being identical to lowering the Mach num- 
ber and/or making the forcing "more solenoidal"). In future work, 
we hope to explore a more detailed approach that includes explicit 
treatment of anisotropy as well as higher-order effects such as am- 
bipolar diffusion; however, this necessitates extensions to our for- 
malism beyond the scope of this paper. 

7 INTERMITTENCY & CORRELATED TURBULENT 
FLUCTUATIONS 

We have thus far assumed that fluctuations on different scales are 
(in the isothermal case) un-correlated and continuous, which leads 
(via the central limit theorem) to locally Gaussian statistics (see 
the discussion in § |2j. However, intermittency implies violations 
of this assumption (manifest in e.g. the deviations of the struc- 
ture functions from self-similar, normal statistics); physically, this 
corresponds to discrete, coherent structures, such as shocks, sound 
waves, and vortices in the turbulent flow. 



0.97 b). Molina et al. (201 2) propose a /^-dependent b correction slightly 
different from that here, which works comparably well. For our purposes, 
both approximations are functionally identical in terms of how they enter 
our formalism. 



However, a wide range of studies have shown that many sta- 
tistical properties of intermittent turbulence can be at least phe- 
nomenologically described by cascade models such as that in She 
|& Levequej (1994) . We briefly review these here. In self-similar 
Kolmogorov turbulence (obeying purely lognormal statistics), the 
velocity structure functions S„(R) = (Sv(R)") = (|v(x) — v (x + 



She 



R)|") should scale as a power-law oc R^" with — n/3. The 
|& Leve que ( 19941 model was originally proposed as an alternative 
scaling law for the structure functions, predicting 



C = (l-7')" + 



1-/3 



n/3 



(35) 



with the original choices j3 = 7' = 2/3 (giving = ;i/9 + 2[l — 
(2/3)" /3 ])HThe meaning of these variables is discussed below 
but this gives a means to parameterize the degree of "non-self- 
similarity" and non-Gaussianity in turbulence. 

|She & Wayimrel <[T995> and |Dubrulle| fl994) showed that 
the scaling in Eq. [35] was the exact result of a general class 
of log-Poisson statistics. They assume extended self-similarity, 
i.e. 5v(R)" oc R" /3 4 /3 (where e R is the dissipation term between 
scales), and a general hierarchical symmetry between scales - such 
that the statistics of e R can be encapsulated in a term n R such that 
e R ~ ttr (where describes the scaling of the average prop- 
erties of the most extreme/singular objects, e.g. one-dimensional 
shocks). Under these conditions, the |She & L eveque ( 1994) scaling 
is equivalent to the statement that ttr obeys log-Poisson statistics of 
the form 



P(n R )d-KR = dY^2P m (m)G R (Y,m), Y = ^ 



with 



Px(m) 



exp (—A) 



(36) 



(37) 



and Gr is any well-defined, infinitely divisible probability distri- 
bution function (physically depending on the driving and charac- 



ter of "structures" in the turbulence). Since 5v(R) oc ej 3 oc 7r 
ki(5v/(<5v)) =ln7r 1/3 = (l/3)ln 7r is a linear transformation and 



1/3 



should obey the same statistics as (l/3)hi7r/; (see e.g. She & 
Way mire 1995, Dubrulle 1994). This describes a general class of 
random multiplicative processes that obey certain basic symmetry 
properties. 

Now recall, the basic assumption behind our approximation 
(and essentially all previous analytic work) for the dispersion-Mach 
number relation and S(R) is that the density field is the product of 
compressive modes in the velocity field, i.e. 8p = p/po — 1 obeys 
the same statistics as 8v(R) (for small "steps" in scale). So, under 
these assumptions, the statistics of In p should have the same form 



as (1/3) ln7r^; and from Eqn. 36 the statistics of (1/3) Itltyr are 
identical to the statistics of Iuttr for a value of /? — > /3 1 ' 3 . Thus, 



»(,W)^&, B (.) 6 (^! 



(38) 



with P P = /3 1 / 3 . We emphasize that this is an assumption: after de- 
riving some basic consequences (§ |7.1[ ( we discuss whether it is 
justified, and how the above formula compares to density PDFs in 
both experiment and numerical simulations (§|7.2|l. 



13 The label 7' owes to historical notation, it is completely unrelated to the 
polytropic index 7. 
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7.1 Monte-Carlo "Steps" In an Intermittent Hierarchy & 
Their Meaning 

Now consider the application of the log-Poisson statistics above 
on the density "'walk" as a function of scale. Following Liu & 
Fang ( 2008 ), first consider the simplest possible form of the driv- 
ing function Gr, a Dirac 8 function (the simplest case considered in 
|She & Waymir^(T995) ; |Dubrulle| ( (T994l l). This corresponds to all 
driving "events," shocks, and other structures which force veloc- 
ity/density changes being strictly quantized and having identical 
fractional magnitudes. In this case, P(lnp) = P\(m = lnp/ln/3 p ). 

In such a log-Poisson random multiplicative process, the (lin- 
ear) variables, in this case pj po, are related by a statistical hierarchy 
between a larger scale Ri and smaller scale R 2 : 



with 



Pr 2 = Wr iRi _ p«, 



WR l R 2 ^f3";(R l /R 2 y 



(39) 



(40) 



Here m is a Poisson random variable with P(rn) = P\ (m) and mean 



A = \r iRi = ln^i/ifa) 
1 p P 



(41) 



from mass conservation. Because the log-Poisson statistics of 
Eq.[36]are infinitely divisible, we can break up the statistics at any 
scale R or \(R) into a sum of "steps" (Ar,r 2 ), and will always re- 
cover the same statistics independent of the step structure. 

Note that this process does not inherently increase or decrease 
the variance in the distribution relative to Gaussian steps. In fact, 
the variable 7' is related by definition to the differential change in 
the variance of the distribution over a "step," b}[^| 



7 



AS 



dS I 1-/3 P 



(42) 



\\a(Ri/R 2 )\ (ln/3 p )2 I din/? I (ln/3 p ) 2 

For clarity, and because the variance-Mach number relation is well- 
established even in systems with large intermittency (see e.g. the 
cases in |Schmidt et aTJ|2008l [20091 |Federrath et~aTl|2008l |20T0l 
|Price & Federrath|2010| > we compare the effects of intermittency at 
fixed variance as a function of scale (i.e. fixed AS). We stress that 
this represents no loss of generality. 

This scaling is sufficient, then, for us to re-construct all of our 
previous predictions, for any choice of /3 P (and 7'), via the Monte- 
Carlo trajectory sampling method. We simply replace the Gaussian 
random step in an interval ln(Ri/R2) with the log-Poisson step 
here, using the same calculation of dS/dlnR and critical density 
p c rit for collapse p^j 

14 This is the same relation (replacing the variance in lnp with that in lnv) 
used to define 7' in the traditional |She & L eveque ( 1994 ) model for the ve- 
locity moments. There, 7' is approximately a constant, because of the ap- 
proximately power-law behavior of v< (R) over the inertial range. Likewise 
here, 7' is nearly constant over the same range. We could, purely mathe- 
matically speaking, hold 7' fixed to some value and either vary /3 to match 
S or allow S to diverge with scale, but this leads to unphysical results. 

15 Strictly speaking, we should also account for intermittency in the turbu- 
lent velocity field, since Vt(R) enters the barrier and variance. However, 
at the level of our approximation only the second moment of v<(ff) en- 
ters, so for a given power spectrum the results are fixed. Accounting for 
higher-order effects would require a multi-dimensional model which can 
account for explicit correlations between the local velocity and density fluc- 
tuations. Numerical si mulations suggest these correlations are small (e.g. 
local v t ~ p~ 04 ; see|Passot & Vazquez-Semadeni||l998| |Kritsuk et al.| 
|2007| jFederrath et al.|201(ty , but they may be important for the most ex- 
treme structures. 



The parameter f3 p = — I represents the degree of intermit- 
tency (and non-Gaussianity) in the distribution: as /3 P decreases 
from 1 to 0, the strength of correlations between fluctuations on dif- 
ferent scales increases and the PDF becomes more non-Gaussian. 
When /3 P — > 1, A — > 00, and the Poisson distribution becomes ap- 
proximately Gaussian and continuous, hence the log-Poisson statis- 
tics become lognormal (this is true for any Gr in Eq. \36\ . When 
P P — > 0, steps are perfectly correlated and the density PDF if in- 
finitely skew. In the structure function derivation of Sh e & Leveque| 
1 1994} , the parameter 7' represents the "degree of singularity" of 
the most intermittent structures on the appropriate scale; as dS/dR 
increases at fixed /3 P these become more singular, raising the vari- 
ance but in concentrated (and locally correlated) structures. 

We should note that the extreme limit where G R is a 5-function 
and the effects of "events" are strictly quantized (as in She & Lev^| 
|eque|199"4) see also |Liu &"F ang 2008) is almost certainly too ex- 
treme - this predicts a quantized (non-continuous) density distribu- 
tion. But it provides a useful upper limit for how strong the effects 
of intermittency might be. In Appendix [B| we consider alternative 
descriptions with continuous Gr, which appear to better describe 
the density statistics in simulations. For our purposes, however, the 
results are similar. 

7.2 Comparison with Simulations and Observations 

We caution that intermittency in the density field of compressible 
turbulence remains poorly understood, and is an important topic 
for future, more detailed study. That said, our assumption mapping 
the intermittent compressive (longitudinal) velocity statistics to the 
density field appears to apply in numerical simulations of super- 
sonic turbu lence ( |Kowal et aijj|2007) |Liu & Fang]|2"008{ |Schmidt| 



etal. 



2009)0 



To test this, in |Hopkins] ( |2012a) , we compile a large suite of 
such numerical simulations (including those above, with and with- 
out magnetic fields, and M ~ 0.1 — 15), and confirm that our "in- 
termittent" density PDF model appears to be a (surprisingly) good 
approximation. We show that in every case considered, a density 
PDF of the form in Eq. [38] (see Appendix [B] Eq. |B4[ l describes 
the simulations accurately (moreso than a pure log-normal), and 
moreover that the "degree of intermittency" (e.g. /3 and 7, or T in 
Appendix [Bj estimated by fitting to the density PDFs is consistent 
with that determined directly from the velocity structure functions. 

This follows from the fact that the density PDF is identical 
to the volume density distribution of a passive scalar in the limit 
of zero diffusion (pure advection) - i.e. mass is a non-diffusive 
passive scalar. The statistics of passive scalars are well-studied, al- 
beit primarily in weakly-compressible turbulence (for a review, see 
Warhaft 2000). Both analytic arguments (from multi-fractal mod- 
els as well as directly from the Navier-Stokes equations; see Yakhot 
|He et al.[l99 8) and numerical simulations (e.g. |Chen & Cao 



1997 



1997 1 have argued for the same conclusion that if the velocity field 



obeys a log-Poisson hierarchy, the PDF of adverted scalars should 



16 |Schmidt et al.| (2009 1 noted large non-lognormal features in such cal- 
culations, and argued that they followed from velocity intermittency, with 
log-density structure functions following directly from the velocity struc- 
ture functions. Kowal et al. (2007) systematically measured the log-density 
PDFs and structure functions in both strongly and weakly-magnetized tur- 
bulence spanning M ~ 0.1 — 7, and found in all cases that the struc- 
ture functions could be well-described by the multi-fractal cascade models 
above. And |Liu IF Fang (20081 showed specifically that a log-Poisson hi- 
erarchy accurately describes the density structure functions and scale-scale 
correlations in simulated cosmological baryonic fluids. 
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as well. This scaling has been measured in laboratory fluid turbu- 
lence for the scalar difference PDF shape vs. scale (the most im- 
portant property here) as well as the scalar structure functions (see 
e.g. |Ruiz-Chavarria et al.fl996l|Zhou & Xia|2010|>. 

Empirically, |Burlaga| ( |1992| l and |Marsch & fu[ ( |1994| l showed 
that the observed PDF, Lagrangian structure functions, and scale- 
correlations of the density fluctuations in the solar wind follow the 
same multi-fractal scaling as the velocity fluctuations (exactly what 
we have assumed here). More recently, Leubner & Voros ( 2 005b|af 
showed specifically that the density statistics agree well with the 
class of log-Poisson models described in Appendix [B] and |Hop-| 
|kins|p012a^. For a r eview, see [C hang (2009}. And in a series of 
papers, |Liu & Fang| p008> ; |Lu et al.| ( |2009F ; |Yuan et al.| ( (2009) ); 
|Lu et al~ 1 2010) have specifically argued that both numerical cos- 
mological simulations and observations of the density field in the 
IGM, cluster gas, and Lya forest appear to exhibit density PDF 
shapes, structure functions, and scale-scale correlations following 
the log-Poisson model. The same model for the density fluctuation 
PDF and structure functions also appears to successfully describe 
laboratory MHD plasma turbulence (seejBudaev 2008 and refer- 
ences therein). 

7.3 Analytic Solutions 

Remarkably, the generalized formulation here still admits analytic 
solutions for the mass function of the nature described in § |2.2| for 
a given assumption about Gr. The full derivation of solutions are 
presented in Appendices |C|D| 



8 RESULTS AT A FIXED "INSTANT" 

Having developed these models, we now examine various "instan- 
taneous" properties of self-gravitating objects in fully-developed 
turbulence (before considering below how these evolve in time in 
§[9}. For the sake of clarity, we will vary several parameters in turn, 
but otherwise fix our "reference model" parameters. This is a rotat- 
ing disk with Toomre 2=1 (marginally stable), large-scale Mach 
number Mi, = 10 (super-sonic turbulence), p — 2 (Burgers turbu- 
lence, appropriate in the supersonic case), b 2 = 1 (i.e. Mi, defined 
in terms of the compressive Mach number) and 7=1, with no in- 
termittency. The dependence on all of these choices is discussed 
below. 

8.1 Basic Properties and Scaling Relations 

First, we consider some basic scaling relations which arise even 
before we evaluate the full density field statistics. 

Fig.|2]plots several such scalings that follow immediately from 
the definition of the critical density and variance S(R). This in- 
cludes the critical density p a n(R) and variance S(R) themselves, 
the mass-size relation for self-gravitating objects, the "barrier" B 
that defines self-gravitating objects at each radius in units of the 
variance {v = fi/S 1 ^ 2 ), and the derivative of the barrier with respect 
to the variance dB/dS. 

In all cases, the variance S(R) on very large scales (^> h) falls 
rapidly; this must be true simply because of mass conservation (fi- 
nite total mass). The variance grows with decreasing R over in- 
termediate scales, where turbulence produces density fluctuations. 
As expected, S(R) increases systematically with Mi, and b (larger 
compressive velocity fluctuations). Lower p (shallower turbulent 
power spectra) contribute more power at small scales, and stiffer 
7 suppress the power on small scales more strongly. Eventually 
the variance converges on small scales below the sonic length (dis- 
cussed below). For non-isothermal cases, the "effective" variance 
S(R, p) is also a function of the density (since the sound-speed, 



hence M is density-dependent) - we therefore plot the variance 
S(R, pent) near the collapse density. For otherwise fixed parameters, 
the variance at the mean density S(R, po) is just the plotted value 
for the isothermal case 7=1. The "turnover" at high-7 corresponds 
to strong suppression of high-density fluctuations on small scales. 

The critical density rises on large scales with support from ro- 
tation («), and on small scales with support from thermal pressure 
and turbulence. Increasing Q linearly increases p a n in turn; p cl -i t also 
increases with a stiffer equation of state (higher 7) owing to thermal 
pressure on small scales. Interestingly, at otherwise fixed parame- 
ters, pcrit is actually lower when the Mach number is increased. This 
is because Q on large scales is held constant, so lower Mi, means 
lower c s /a g (h), i.e. a lower threshold where thermal pressure be- 
comes important. 

These scalings highly two critical scales in this problem. 

8.1.1 The Maximal Instability ( Quasi-Toomre ) Scale 

On large scales, there is a clear "most unstable scale" 7? ma ximai in 
Fig. [2] This is closely related to the minimum in p cr j t at R w h, 
with rotation making the system more stable on larger scales, and 
thermal support "stabilizing" smaller scales (as in the traditional 
Toomre stability analysis). However the "most unstable scale" is 
not simply where p cr u is minimized, nor is it the Toomre wavelength 
(most unstable wavelength in a constant-p disk). 

To estimate the scale which will contain most of the collapsing 
mass in a non-uniform density disk, the more relevant parameter is 
the ratio v = B/S [ ^ 2 . If the density PDF is lognormal, then 



F co n = -Erfc 



2 LV2J 



(43) 



is the fraction of mass which is self-gravitating on that scale (in- 
dependent of the exact size/mass spectrum), so F co u is maximized 
where v is minimized. We see in Fig.[2]that the rise in p cl -i t on small 
scales is sufficiently rapid that the scale where v is minimized is 
always close to R ~ h (with weak dependence on any parameters). 

Although there is no trivial analytic solution, it is easy (if de- 
sired) to numerically solve for /? ma ximai = R[v mm ]. And if we assume 
turbulence dominates over thermal pressure at the scale of interesf^] 
and expand about R[i/ m m] ~ h, we obtain 



, / ^maximal \ 



A, = 



Ai - y/A 2 + 2A 2 In (2Q [1 +gj) 
2A 2 

2(3 -p) + (1 + k 2 ) (|dS/dlntf|,, - 1) 



(44) 



3 3k 2 (3- P ) 2 
2 16 4(1 + k 2 ) 

This scaling almost always gives R[v mm ] ~ h, but includes the 
weak dependence on other parameters. Higher Q values and steeper 
power spectra systematically shift i? max imai to smaller scales (they 
raise p cr it(^ = h), so require going to smaller R, where S(R) is 
larger, to reach the minimum in v). 

These effects, while weak, are not captured in the traditional 
Toomre analysis. In fact a Toomre-style linear stability analysis of 
a thin disk with constant density gives the most unstable wave- 
length R/h ~ |feni n «| _ = (using h = <j/Q). The Toomre 
wavelength at the mean density increases linearly with Q; how- 
ever, the "most unstable" wavelength - in terms of where most of 



17 If thermal support dominates near ~ h (e.g. Ait, <C 1), we derive an 
expression for Rmaximal which is identical to Eq. [44] with p = 1. 
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Figure 2. Some basic properties of self-gravitating regions in turbulent media. Each assumes a disk with (arbitrary) mean midplane density po, scale-height h, 
and the following "default" properties, each varied in turn. (1) Large-scale Mach number A^/, = Mtm{h) = 10, varied from 0.1 — 30. (2) Turbulent spectral 
index (E(k) oc k~ p ), p = 2 for Burgers (highly super-sonic) turbulence is our defaul t an d p = 5/3 is Kolmogorov, values outside this range are rare. (3) 
Global Toomre parameter (default Q = \, for marginal stability). (4) Parameter b in Eq.[l2](the Mach number-density dispersion relation), b 2 = 1 (default) for 
Ml, defined as only the compressive part of the turbulence; more generally b 2 = 1/3 in strongly magnetized media, b 2 = 1/4 in randomly forced supersonic 
turbulence, b 2 = 1/9 for purely solenoidal forced turbulence. (5) Equation of state polytropic index 7 (isothermal gas with 7 = 1 by default, which produces 
a lognormal density PDF; other choices lead to non-lognormal distributions). Rows show different properties. Top: Mass-radius relation; this follows Eq. |47| 
and generically produces weakly-varying surface densities over a wide dynamic range. Second: Variance S(R) in the density field; for non-isothermal cases 
(7 7^ 1) we plot S(R, Pcrit)> the variance at R for a narrow density range about the critical density. Larger M/, (and b) systematically increase 5, stiffer 7 
suppresses it on small scales. Third: Critical density (in units of po) f° r self-gravitating collapse. This increases with Q and 7 (on small scales); but at fixed Q 
it is lower at higher Mh- Fourth: dB cr i t /dS (B C rit is the collapse threshold or "barrier" in Fig.[TJ; when 2> 1, the collapse threshold rises "too fast" relative to 
density fluctuations and the field becomes non-self gravitating (the "last-crossing regime"), when > but <g 1 regions are self-gravitating up to large scales 
("first crossing" regime); when <C on large-scales, direct collapse is suppressed and structure growth switches from "top-down" to "bottom-up" via mergers. 
Bottom: u = B CIlt / 'S 1 / 2 ; the number of standard deviations of a collapsing fluctuation on that scale. Minima here correspond to the most unstable scale (~ h). 



the mass will fragment under self-gravity - actually decreases with 
increasing Q in a turbulent medium. 

8.1.2 The Sonic Scale 

Another critical scale is the sonic length/mass scale. For isother- 
mal gas this is typically defined as the scale below which the rms 
A4 < 1, Rf, mc ~ R(A4 — 1). The corresponding "critical mass" 



is M S onic = M(.K son i c ). But recall, what matters for the density dis- 
tribution is really the compressive component of the Mach num- 
ber .Mcompressive, which enters the variance as M compressive = bM, 
hence the effective sonic scale of interest for the mass function is 
R(bM = 1). 

For the polytropic case, there is no one spatial scale where 
bM. = 1, but since we care about objects crossing the critical 
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Figure 3. Mass functions (MFs) of self-gravitating objects, for the same parameter variations as Fig. [2] We plot Mdn/d log M, the mass per logarithmic interval 
in object mass (a MF with dn/dM oc M~ a with a = 2 is a horizontal line). Thick lines show the first-crossing distribution (mass spectra of self-gravitating 
regions defined on the largest self-gravitating scale). These characteristically have faint-end slopes a < 2 and the mass is concentrated near the "maximal 
instability" scale (collapse on scale ~ h), with weak dependence on other properties except at low masses. Thin lines show the last-crossing distribution (mass 
spectra on the smallest self-gravitating, i.e. non-fragmenting, scale). These have slopes slightly steeper than a > 2 over a broad range, then turn over below 
the sonic scale, defined by Eq. |46| and labeled in each panel by the short vertical line (when Mi, < 1, this is not defined, but the relevant scale is now again 
the maximal instability scale). The normalization of the MF decreases for Mi, < 1 or i) < 1 but is otherwise similar (saturated when an order-unity mass 
fraction is self-gravitating). The dynamic range of first & last crossings increases with M max j ma i/M SOI1 i c . This is a strong function of Mi,, modest function of 
Q and p (though over the expected range p = 5/3 — 2 the effect is small) and weak function of b (except i)< 1). Non isothermal 7 have no effect in the 
turbulence-dominated scales > M son j c , but increasing 7 sharply suppresses first and last crossings below M son j c , with the cutoff becoming steeper with higher 
7 (infinitely steep for 7 > 4/3). Intermittency has small effects on the MFs (see Appendix|A|B|. 



density, we should consider M at this density, giving /?sonic = 
R(b _M[pqrit(2?)] = 1)- Although the general solution for p a i t and 
Rsomc for arbitrary 7 must be numerical, it is possible to analyti- 
cally solve for i? son i c and Af SO nic if we assume either R son i c <C h or 
non-rotating turbulence (i.e. neglect the large-scale/«: terms, gener- 
ally a good approximation if Mi, > !)■ This gives 



equations of state it implies more turbulence on small scales, and 
lower _R sonic . 

Below iconic, the contributions to the density variance from 
successive turbulent fluctuations are small. This scale is obvious in 
Fig. [5] as the scale below which S(R) saturates and /9 cr j t begins to 
steeply rise, producing a very sharp increase in AB/AS. 



2Q' j Li 1 



M SI 



Ma 



{bMci 



(bMci, 



m = 



V Ra 

7-1 



3-H 



(45) 
(46) 



2(7-l)-(p~l)(7-2) 



For Aid S> 1 and p = 2, the case of greatest interest in the ISM, 
this reduces to (R s<mic /R c {) w (22') (7 "' )/7 fc (27 " 4)/7 M^ h and 
(M somc /M cl ) « (20') (37 - 2)/7 6 (67 - 8)/7 M7 4/7 . 

Stronger turbulence (higher Mi,), softer equations of state 7, 
and shallower (more bottom-heavy) turbulent power spectra (lower 
p) decrease R SO m C and M son j c as expected. The sonic scale is inde- 
pendent of Q for isothermal gas but for "stiff" equations of state, 
higher Q implies more support and higher ^ son i c while for "soft" 



8.1.3 Mass-Radius Relation, Densities, & Surface Densities 

From p a -,t(R), it is trivial to calculate the corresponding mass-size 
relation of self-gravitating objects at their formation (Eq. |16} . We 
see in Fig. [2] that these obey approximate power-laws over a very 
wide dynamic range, albeit with some gradual shifts in slope. 

The critical densities p cr i t shown in Fig.|2]rise at large scales 
because of angular momentum support (k), then rise again below h 
from a combination of turbulent support (above iconic) and thermal 
support (below) which must be overcome by gravity. Over each 
regime where thermal, turbulent, or rotation support dominates, we 
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AM = log(M last /M first ) AM = log(M last /M first ) AM = log( M last / M firsi ) AM = log( M last / M first ) AM = log( M last / M first ) 



Figure 4. Dynamic range of the "fragmentation cascades" that make up the MFs in Fig. [5] Top: Probability distribution, for a given random point in space 
(i.e. random volume element) which is self-gravitating on some scale, of the logarithmic dynamic range between the scale of first crossing i?g lst and scale 
of last-crossing i?i ast . Most volume within a self-gravitating region is only self-gravitating over a narrow dynamic range, i.e. is not part of a much denser 
fragmenting substructure. The "tail" extends in dynamic range with the ratio M somc /M max j ma i, as the gap between first/last crossing in Fig. [3] Bottom: Same, 
but for a random self-gravitating mass element, as a function of the gap in first/last crossing mass scales. This is distributed over a wider dynamic range (note 
the difference in figure scales). A given Lagrangian element is often self-gravitating over several dex in mass scale (i.e. locally self-gravitating structures prefer 
much larger "parent" structures). The predicted dynamic range generally scales with M somc /M max j ma |. However varying p produces weaker (or even opposite) 
effects than what might be expected: when p is "flatter" (more turbulent power on all scales), the last-crossing MF and M sonlc do extend to smaller masses, but 
these often correspond to smaller-scale first-crossings as well (rather than being entirely seeded inside much larger structures). 



can determine the approximate power-law: 



8.2.1 First Crossings - The Largest Collapsing Objects 



M(R) 
po h 3 



R\3 



(R > h) 
(47) 



In the turbulence-dominated regime (^sonic R h\ the size- 
mass relation is simply set by the turbulent power spectrum. For 
typical turbulent properties, p ~ 5/3 — 2, this generically gives 
nearly-constant surface densities E ~ M/R 2 in collapsing objects. 
In the thermal-dominated regime (R <C ^ SO ni c ), the fact that M(R) 
is multivalued for 7 > 4/3 is intimately related to the suppression 
of fragmentation: if an object attempted to fragment or contract to 
smaller spatial scales, it would have to acquire more mass to over- 
whelm the rapidly increasing thermal pressure barrier. 



8.2 The Mass Function: Gravitating Structures and 
Sub-Structures 

Fig. [3] shows the mass functions of both "first-crossing" objects 
(self-gravitating objects with masses/sizes defined on the largest 
regions on which they are self-gravitating) and "last-crossing" ob- 
jects (defined on the smallest self-gravitating scale). 



The first-crossing MF has most of the mass concentrated near the 
maximal instability (quasi-Toomre) scale defined above (the scale 
most vulnerable to fragmentation). This is close to ~ h in all 
cases, but with a weak dependence on other parameters (decreas- 
ing weakly with Q, the opposite of the behavior of the traditional 
Toomre length/mass; see § |8.1.1[ >. 

On larger scales, angular momentum produces an exponential 
cutoff of the MF (as the variance falls of rapidly and v grows with 
R). On smaller scales, there is a power-law like run over ~ 1 — 4 
dex in mass, the range where turbulence dominates. As discussed 
in Paper I, if we use the analytic first-crossing solution for the linear 
barrier to approximate this, we can indeed predict the generic slope 
An/ AM oc M~ Q with a slightly smaller than 2: 



"first < 



(3-P) 2 
2p 2 S 



In 



M 



aximal 



(48) 



The value a ~ 2 is generic because gravity is scale-free, so 
if the critical density and variance per unit scale were constant we 
would recover exactly equal mass per logarithmic interval in mass 
(a = 2). These quantities are not, however, constant; but they are 
close to it: B and S are logarithmic functions of R and M, because 
the density distribution follows a random multiplicative process. 
As a result, there is only a weak correction: because v increases 
towards smaller R (see Fig. [2} the MF has less mass at small-M 
than a = 2 (i.e. we predict a < 2). 

At sufficiently small scales, the MF again cuts off. We expect 
this to occur as we approach the sonic mass, since density fluctua- 
tions are strongly suppressed (Fig. [TJ; plotting this in each case in 
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Figure 5. Correlation function £(r) (excess probability of a neighboring 
object within a distance r of a given object, Eqs. |21|22) for last-crossing 
(left) and first-crossing (right) objects of a given mass. Left: Last-crossings 
with mass M[ ast (defined as the mass at the maximum of the MF; the peak of 
Mdni as t/dlnM in Fig.[5J. In general Mi ast ~ M son i c . Note £(r) is undefined 
below r pb R(M\. ist ). Right: Same for first-crossings at Mfl rst (the peak of 
Md«fl rs t/dlnM), roughly Mfi lst ~ M max ; ma i. For each, we consider a range 
of parameters. Measured at these "characteristic" masses of the MF, £(r) is 
nearly universal. There is a dependence on M/, when Mi, < 1 ; here £ (r \ M) 
increases as the number density n(M) decreases. The shape is generically 
power law-like, with §(r) <x r~ 2 at large separations transitioning to shal- 
lower oc r~ 1 at smaller separations. 



Fig. [3] we see this is indeed a good predictor of the low-mass MF 
cutoff. In all cases, masses near M son j c contain only a small fraction 
of the total mass in first-crossing objects. 

Interestingly, the first-crossing MF slope (in the power-law 
range) is remarkably independent of the turbulent spectral index p, 
despite the fact that this sets the run in M(R) and p a i t (R)- It appears 
that the steeper run in 5 (steeper at higher p) nearly completely can- 
cels the steeper run in p, giving the quite similar behavior in v(R) 
in Fig. [2] as well as a nearly identical dfi/dS from ~ h down to 
Rsamc' This is not trivially expected from the first-crossing solution 
for a linear barrier, and depends on the full (non-linear) solution 
of the MF. The index p does make some difference at low masses, 
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Figure 6. The correlation function §(r) as Fig. [5] but for objects with fixed 
mass M. Left: £(r) for last-crossings with mass M = 0.01 pa /1 3 . Right: £(r) 
for first-crossings with mass M = pair' ■ The dependence on other parame- 
ters is still weak, but stronger than in Fig. [5] This dependence comes from 
sampling different parts of the MF relative to the peak (sonic or maximal 
instability) mass scale. The amplitude of £(r) systematically increases with 
decreasing number density at fixed M (Eq. |51| . The last-crossing £(r) con- 
verges to universal shape at large r, but becomes more shallow at small 
radius when the variance is larger (explaining the trends in Q, h, p, and 7). 



where it defines how rapidly the transition from turbulent to ther- 
mal support occurs. 

The MF does depend significantly on Mi, and b. Increasing 
these increases the variance at all scales and so broadens the range 
where collapse is predicted. The effect is weak near the peak of the 
MF in the supersonic case and/or when b 2 > 1/4. In these cases the 
peak of the MF is "saturated" (i.e. reaches near order-unity frac- 
tions of the mass), so its normalization changes negligibly even 
if we vary the properties of the turbulence dramatically; moreover 
the power-law regime is in the near-universal range so it also does 
not vary much. But increasing Mi, has a large effect on iconic and 
Msonic, correspondingly extending the MF to lower masses with in- 
creasing Mi,. For Mi, <C 1 or b <C 1, the MF begins to be sup- 
pressed rapidly. In this limit, the variance in the density field is 
small, so v becomes 2> 1 even near the most unstable scale. As 
such, we transition from sampling the "core" of the density distri- 
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Figure 7. Global, time-averaged rate of formation of bound structures, as a 
function of mass, for a model in which all density modes start at 8 = (me- 
dian density) and self-gravitating structures are immediately removed (i.e. 
collapse or are destroyed). We show both first (black) and last (red) cross- 
ings. Histograms show the full numerical result (for the "default" model). 
Dotted lines show the instantaneous MF predicted for a fully-developed 
density PDF (Fig.[5J, divided by a constant timescale fcross (.Rj?^) (the tur- 
bulent crossing time at the scale where the MF peaks in Fig. [5). Dashed lines 
show the same "fully-developed" MF, divided by the crossing time at each 
mass/radius scale f C ross(^) (for first-crossing we also offset it by 0.3 dex 
to lower masses). This provides a good approximation to first-crossings: 
they develop approximately independently on each scale on the local cross- 
ing time, but with a truncated high-mass MF (because the highest masses 
depend on merger-driven growth, discussed below). However the constant 
timescale (dotted) is a better approximation for last-crossing: these are not 
driven by rapidly-evolving, small-scale fluctuations, but "pre-seeded" in- 
side larger objects and so tied to the rate of first-crossing events driven by 
larger-scale fluctuations. 



bution to the tails, and the probability of collapse is exponentially 
suppressed. Because, within the tail, the probability of a "cross- 
ing event" is a strong function of v, the MF also rapidly becomes 
sharply peaked around the most unstable scale. 

The MF also depends on Q. For Q < 1 , the fact that the system 
is unstable even at p — po means that essentially all of the mass is 
in self-gravitating objects on the maximal instability scale (making 
a much sharper peak). For 2 > 1, we again move into the expo- 
nentially suppressed regime as with Mi, <C 1. 

Unsurprisingly, over the turbulence-dominated regime, there 
is essentially no dependence of the MF on 7. Below the sonic mass, 
though, 7 makes a large difference, with more mass in low-mass 
objects with softer 7. 

8.2.2 Last Crossings - The Smallest Collapsing Objects 

The last-crossing MF is also truncated exponentially at high 
masses; this always occurs at somewhat smaller masses than the 
peak of the first-crossing MF. There is also a near-exponential cut- 
off below the sonic mass. 

In general, the last-crossing MF is "saturated" when Mh > 1, 
and the effect of changing various parameters is largely to shift the 
sonic mass/lower cutoff of the MF. However as for first-crossings, 
large Q 3> 1, Mi, <C 1 or b <C 1 lead to an exponentially suppressed 
MF with a more narrow mass range around the most unstable mass. 
This forces convergence between the first and last-crossing MFs. 

There are some generic differences from the first-crossing MF. 
There, the MF generically had a slope shallower than a = 2 (i.e. 



had most mass in the most massive objects); here the slope is much 
closer to a ~ 2 and even slightly steeper over some range. 

Mathematica lly, th is is related to the exact solutions for a lin- 
ear barrier in Eq. 4 18 The only difference between first and last- 
crossing MFs in that case is the replacement of the term Bo/S in 
front of the first-crossing MF with dfi/dS for the last-crossing MF. 
Bo and dB/dS are constant by definition for a linear barrier, but 
5 = S(R) decreases with R, so the first-crossing MF should be shal- 
lower than the last-crossing MF by a term oc S~ [ . From Fig. [2] we 
see this is a weak logarithmic power ~ M°' 2 ~ 0A . 

Physically, the origin of this difference is that, when there is 
close to equal power on all scales, the collapsed fraction should 
be close to uniform on all scales; but the first-crossing MF only 
"counts" objects on the largest scales in which they are self- 
gravitating. Most of the self-gravitating mass on small scales lies 
within larger objects, rather than being isolated. This is closely re- 
lated to the correlation functions discussed below. 

We see a stronger dependence on 7 in the last-crossing MF 
as compared to the first-crossing MF, because more of the mass is 
concentrated near the sonic mass (where the gas thermal physics 
is important). For 7 < 1, the low-mass cutoff is significantly more 
shallow; for 7 > 1, it becomes quite sharp. For 7 > 4/3, collapse 
is entirely suppressed below /f S omc> so there is a sharp "spike" then 
cutoff in the MF. 

In the last-crossing MF we see a stronger dependence of the 
MF slope over the turbulence-dominated range on p. The depen- 
dence is still weak; but as shown in Paper II, if we use the linear 
barrier exact solution and drop the angular momentum and thermal 
pressure terms, we expect a slope 



«liist c 



(3 -pf In (M/M m 



- p In 2 



2S(M)p 2 



(49) 



for the values of M and 5 here, this becomes slightly shallower than 
a = 2 at lower p < 2, and steeper at p > 2. However over the most 
likely range p ~ 5/3 — 2, the difference in slope is small (~ 0. 15). 

8.3 The Dynamic Range of Fragmentation 

From the comparison of the MFs in Fig. [3] it is tempting to in- 
fer the "dynamic range" of fragmentation, i.e. the range of scales 
over which fragmentation occurs. This should be crudely related to 
the "gap" between the peak of the first-crossing MF and the sonic 
scale cutoff in the last-crossing MF. But the MFs alone to not tell 
us how an individual "trajectory" (Eulerian or Lagrangian volume) 
behaves, in particular the range of scales over which it (individu- 
ally) is gravitationally unstable and can fragment. 

We calculate this in Fig. [4] To do so, we use the Monte Carlo 
method described in § [2] and for each parameter set from Fig. [3] 
generate ~ 10 s "trajectories" sampling the volume. For each trajec- 
tory that is anywhere self-gravitating, we then record the location 
of first and last crossing, and measure the logarithmic interval in ra- 
dius AR = log (i?iast/^flrst) over which that volume element is self- 
gravitating. Since each trajectory represents a volume element, this 
corresponds to the distribution of spatial scale ranges over which 
fragmentation will occur, for a given random volume (i.e. random 
point in space) selected in any self-gravitating object. 

We can also compare the distribution of mass scales over 
which fragmentation will occur: we follow the same procedure, but 
now define the "mass range" AM = log (M\ ast [R\Mt]/Mn Ist [Ra Ist ]), 



18 We can see in Fig. [2] that dB/dS ~ constant over much of the dynamic 
range of interest, so this is not a bad approximation to the full solutions. 
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Figure 8. The steady-state MF which develops in time from uniform (S = 0) 
initial conditions. The density field evolves and crossings are recorded as 
Fig.|7] but regions which cross (become self-gravitating) are allowed to sur- 
vive and evolve for some finite time fbound °c kmss (R) (with R evaluated at 
the time of first crossing). For fbound <C fcross, the steady-state MF resem- 
bles the instantaneous rate-of-formation (more sharply concentrated near 
the most unstable mass) multiplied by a lifetime (bound- As (bound increases, 
structure develops "top-down" at low masses (building up the low-mass end 
of the MF) and "bottom-up" at high masses. As fbound — > oo, we recover the 
analytic prediction for the instantaneous MF in fully-developed turbulence 
(dotted lines). 



and mass-weight so that the distribution reflects the PDF for a ran- 
dom mass element. 

The dynamic range in mass and spatial scale of fragmenta- 
tion behave broadly as expected: they become more broad with 
increasing Mi„ lower Q, and softer 7, all in line with the lower 
sonic lengths/masses predicted (while the upper mass, at the max- 
imal instability scale, stays about constant). Larger b (more highly 
compressive turbulence) also increases the dynamic range, as seen 
in the width of the MFs, by extending the range over which the 
variance is significant compared to the barrier. 

Interestingly, even though the MF is more broad with lower p 
(shallower turbulent power spectra), the "typical" dynamic range of 
fragmentation is actually smaller. Shallower power spectra spread 
the turbulent power over a wider range of scales more uniformly 
(with p = 1 corresponding to equal power on all scales), making 
it more likely that a fluctuation can become self-gravitating (cross 
above and below the barrier) on small scales, without necessarily 
being embedded in a more powerful larger-scale fluctuation. Thus 
even though the MF is more broad, the local dynamic range of frag- 
mentation within each self-gravitating region is not. 

In a mass-weighted sense, the distribution of fragmentation 
scales is relatively flat, over the range AM ~ log (M son ic/M ra aximai), 
corresponding to the scale-free MFs. But in the volume-weighted 
sense, the distributions are significantly more narrow. Most of the 
volume, selected within a self-gravitating object, is not likely to 
be self-gravitating on much smaller scales. Fragmentation occurs 
in dense sub-regions that occupy a small fraction of the "parent" 
volume. 

8.4 Correlation Functions 



Figs. |5|6| show the predicted three-dimensional autocorrelation 
functions (f (r) from Eqs. 21|22 1 for first and last-crossing objects, 
as a function of the parameters varied in Figs. |2|3| above. 

We first consider £(r) for objects with the "characteristic" 
mass at which the mass density Md;?/dlogM in the correspond- 



M ol (t=0) = 0.1 



o 
c 

TD 



o 



t=1 6.1 
t=16.7 
t=17.0 
t=17.2 
t=1 7.3 
t=17.3 
t=17.4 
t=17.4 
t=17.4 
t=17.4 



<M C 

u c h)= 

3vl 



0.6 
0.8 
1.1 
1.5 

2.0 

4.9 
6.7 
9.0 



1, 



-3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 



o 
c 

TD 



O 




-2.0 -1.5 -1.0 
log( M / M d ) 



0.0 



Figure 9. Fragmentation mass spectrum (i.e. last-crossing mass spectrum) 
developed as a function of time inside an initially non-fragmenting, self- 
gravitating cloud (i.e. within a previous last-crossing "region") as that re- 
gion collapses (following § |9.2| at constant virial parameter). We show the 
MF at several times (each labeled f in units of ?o (Eq. [63), with the cloud- 
averaged Mach number A4 c {(t) increasing as the region collapses), until 
the entire volume has become self-gravitating on some sub-scale. Mass is 
in units of the total cloud mass M Q \. We compare our standard (isother- 
mal) model, but consider a cloud with initial cloud-scale Mach number 
M c \(t = 0) = 0.1 (top) and 10 (bottom). The former undergoes essentially 
no fragmentation until a very large number of dynamical times have passed 
and the cloud has contracted to where M c \ > 1 — 3 (collapse by a fac- 
tor R c i ~ 0.1 Rq). Almost all fragmentations occur near the "parent" cloud 
scale. A cloud which begins highly super-sonic (despite having no initial 
fragmenting sub-regions) develops fragmentation almost immediately as it 
collapses. 



ing last/first-crossing MF peaks. This is generally close to M son j c 
and M m a X imai for last/first crossings, respectively. With this choice, 
the correlation functions are nearly universal. Their normalization 
depends only weakly on varied parameters, and the shape is almost 
constant; £(r) is a near power-law, but with some run in slope from 
a shallower £ (r) oc at small separations to a steeper £(r) oc r~ 2 
at large separations. 

If we instead consider a constant mass at which to evaluate 
£(r|Af), we see a stronger dependence on the varied parameters; 
the change in clustering reflects the "position" in the MF relative 
to its peak, then, rather than an intrinsic difference in the clustering 
structure around the characteristic mass. 

We can gain considerable insight from the closed-form solu- 
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tion for a linear barrier, derived in Paper III: 



&M(r,M\B = B + nS)- 



exp 



(b 2 m IS, 



M + S 2 M /S(r)]- 



y/l - [S(r)/S M ] 2 



1 (50) 



where S M = S(M) = S(R[M]) and B M = B(M) for a last-crossing 
of scale M[R]. Comparing this to the solutions in Figs. |5|6| shows it 
is a very good approximation. 

At large separations r (where S(r) <C 1), this becomes 
£(r|M) ~ (Bm/Sm) 2 S(r) = u M S(r). Since we calculate £ at the 
characteristic mass of the MF, we are implicitly choosing a point 
near the minimum in v, so vu ~ v mm ~ 1 in each case, hence we 
obtain a similar normalization for But when we choose dif- 
ferent masses, we systematically alter the normalization with vm\ 
more rare systems (higher-^) give a correlation function with sys- 
tematically higher normalization. In this large-separation limit, the 
shape of £(r) simply traces the variance as a function of scale in 
the density field, with some normalization "bias" (that is identical 
in meaning here to the bias parameter of cosmological correlation 
functions). At large scales (above the injection scale or characteris- 
tic scale h set by the maximum velocity dispersions) this declines 
as r~ 2 because S(r) is suppressed by the angular momentum bar- 
rier in Eq. [12] (Fig. |2j a falloff of this nature if required to ensure 
mass conservation). 

At separations r <C h, Eq. |50] becomes ~ 
exp (u M /2)/ ^2 (1 — S(j) /Sm)- The normalization again is 
just a function of um- In Paper III (Eq. 34 therein) we show 
that expanding this about the sonic scale gives an approximate 
power-law £(r) oc r x over the turbulence-dominated intermediate 
regime, where \ smoothly varies from —2 at large separations to 
a smaller value ~ — (13 + p (p — 6))/ (4(1— p)) below ~ h. So 
the slope of at r < h is set by the turbulent power spectrum 
controlling the run in the density fluctuations with scale, and 
is systematically steeper at lower-/?. Finally £(r) steepens as it 
approaches the sonic length and S(r) — > Sm, because there is little 
range for independent random fluctuations to act. 

In Fig. [5] we only see a large change in when the global 
Mach number A4i, is subsonic. This is also where we saw the MF 
in Fig. [3] become exponentially suppressed. This relates to the vm 
dependence discussed above. Recall, for the same linear barrier that 
predicts Eq.[50] we have fi(M) = (fj,/\/2n Sm) exp (— v M /2). So at 
scales r < h (where S(r) ~ 5m), we expect £(r | M) oc exp (y M /2) oc 
fe(M)~\ i.e. the clustering amplitude is inversely proportional to 
the predicted number density! 

Indeed most of the mass dependence and differences in 
in Fig.|6]can be explained by a near-universal number density-bias 
relationship, specifically (for a linear barrier) 



3m(r\M): 



^2Tv[SM-S(r)] 
|dS M /dlnM| 



V M 



dn(M) \ 
dlnM / 



(51) 



V2vr[S M -S(r)] 

with V M =M/p ait {M). 

Physically, fragmentation on small scales is highly clustered 
because most of the power in the density fluctuations comes from 
large scales (see Paper III). Near the sonic length (the characteris- 
tic scale of the last-crossing distribution), the probability of an in- 
dependent fluctuation suddenly arising, localized on small scales, 
sufficient to cross the barrier, is extremely small. Indeed the fact 
that it is the characteristic "last-crossing" mass implies that such 
fluctuations must be vanishingly rare. Thus last-crossings depend 
on fluctuations on larger scales to "seed" most of the necessary 
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Figure 10. Distribution (mass-weighted) of instantaneous Mach numbers 
jVlcl at which an initially non-fragmenting, isothermal self-gravitating 
cloud will fragment as it contracts (as in Fig. [9). We consider the two cases 
in Fig.[5]and plot the average fraction of mass which would fragment and 
collapse per logarithmic interval in M a \ (which is directly related to the 
time since the cloud began to contract and/or total contraction factor), as 
well as two additional cases. Clouds which begin sub-sonic all contract un- 
til Alcl ~ 1 — 3 (independent of the initial M c {) then fragment. Clouds 
which begin with A4 C \ 2> 1 begin to fragment very rapidly. 



density fluctuation, i.e. they "live" inside larger collapsing objects 
(as with e.g. protostellar cores seeded inside of GMCs). This is 
true for any plausible bulk properties of the turbulence (any turbu- 
lent spectral index p > 1). The probability of multiple sub-objects 
forming inside a given region is largely set by the amplitude of the 
larger-scale fluctuation on the scale of that region. 

For this reason, the shape of the correlation function, as we 



see in Eq. 51 is almost entirely set by the run in S(r), which in the 



proper dimensionless units (rescaled between the maximal insta- 
bility scale and sonic scale) is nearly universal (with only a weak 
logarithmic dependence on the shape of the turbulent power spec- 
trum or polytropic index). 

The amplitude of clustering depends inversely on the average 
number density, or more specifically from Eq.[5TJfhe volume filling 
factor, of objects of a given mass, since the higher this is, the "eas- 
ier" it is for totally independent fluctuations localized on a given 
scale to produce collapsing objects of a given mass (hence it is less 
dependent on parent-scale fluctuations). 

8.5 Non-Gaussianity, Intermittency, and Correlated 
Turbulent Fluctuations 

In Appendix[A] we consider in detail how different models for non- 
Gaussian and/or correlated structures in the density field affect the 
mass and correlation functions predicted here. This is parameter- 
ized in our models via either our approximate intermittency treat- 
ment (§|7j or the non-Gaussian and correlated fluctuation structure 
when 7 7^ 1 (§[3}. 

For realistic intermittency levels and equations-of-state, we 
find that the effects on the mass functions stemming specifically 
from non-Gaussian statistics and complicated correlations in the 
density field are relatively minor (generally equivalent to small 
changes in the Mach number or Q; see Figs. |A1|A4[ |. We see in 
Figs. |5|6| that the correlation functions predicted are also not very 
sensitive to the inherent correlation structure of the different scale- 
modes in the turbulence. 

For a detailed discussion, see § [A] But this is generally true 
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because the most important parameters in our model (the width of 
the density distribution, and the run of this width with scale), fol- 
low from the lowest-order moments of the density field - they are 
essentially set by the power spectrum, which varies only weakly 
with these higher-order effects. 

We stress that even our simplest models (7=1 and j3 = 1), 
where we assume uncorrected Fourier fluctuations, produce large 
non-zero correlation functions. In other words, the assumption of 
uncorrelated fluctuation phases in Fourier space is not equivalent 
to an assumption of no real-space correlations in the density field. 
In fact, if we directly calculate the mass or density autocorrela- 
tion function from these models, we obtain a prediction in good 
agreement with that directly measured in full numerical simula- 
tions in Vazquez-Semadeni & Garcia (2001), without needing to 
invoke any additional correlation structure in the field. The same 
is true for the well-studied cosmological case: models with un- 
correlated Fourier modes still produce a real-space mass auto- 
correlation function with significant structure! 

9 TIME DEPENDENCE 

Thus far, all of our calculations have concerned instantaneous prop- 
erties at a fixed moment in time (assuming fully-developed turbu- 
lence). In this section, we extend the models to approximate the 
time-dependent evolution of turbulent density fields and fragmen- 
tation cascades. 

9.1 In a Stationary Background 

First consider a system with constant mean background conditions 
(average density and steady-state turbulent power spectrum). 

If the density field is strictly lognormal with un-correlated 
fluctuations on different (Fourier) scales (as we assume in our sim- 
plest model), then we show in Paper I that the evolution of the 
modes of the density field must obey a generalized Fokker-Planck 
equation. For a (zero mean) Gaussian variable x with variance S x 
the PDF to find the sv_stem with value x at time t given an initial xo 
atfo = t — At is thei 



s yst 



P(x, t) Ax — 



1 



= exp 



:-xo) 2 
2S X 



dv 



(52) 



& = S,[l-exp(-2p-fo]/r)] 
x = x exp(-|> — to]/r) 

Here, r is the correlation tim^J equivalently, it is the 
timescale at which the variance in x(t) with respect to xq grows, or 
the timescale for turbulent mixing of a non-diffusive passive scalar 
(see § \7.2\. This is measured in numerical simulations (see e.g. 



|Klessen||2000l |Ostriker et al.|2001| |Kitsionas et 1012009) , which 
show it is approximately the crossing time 



= ¥R/(v*(R)) 1/2 



(53) 



with f « 1 a constant \Pan & Scannapieco|2010|fin d f « 0.90 — 
1 .05 over the range M ~ 1 .2 - 10; |Zhou & Xia|20l0l obtain similar 
results from experiment in sub-sonic cases). 



" As shown in Paper I, in the limit of small At, Eq. [52] represents the 
only form that the evolution of P(x, t) can take if we require that P is ex- 
actly Gaussian, modes are uncorrelated, S x is conserved in ensemble av- 
erage (true if the turbulence is steady-state), and the growth in variance 
between x and xq is independent of our choice of integration step-size. 
20 By definition, the Lagrangian correlation amplitude between modes de- 
cays with e-folding time r. 




1 10 
Mach Number M d 

Figure 11. Distribution of Mach numbers at which an initially non- 
fragmenting self-gravitating cloud will fragment, as in Fig. |10| for other- 
wise identical initial conditions but varied model parameters as labeled. All 
clouds here begin with M z \ = 0.3, but all sub-sonic initial conditions give 
nearly identical results for each parameter value (and initial _M c i 2> 1 al- 
ways produces very rapid fragmentation, as in Fig. |10( . Changing the power 
spectrum p has weak effects, though very shallow spectra show slightly 
slower fragmentation. Changing the value b of compressive-to-total fluctu- 
ations has the largest effect; collapse requires supersonic compressive Mach 
numbers, i.e. bM c { > 1 — 3. Changing the polytropic index to either softer 
or stiffer values leads to a more narrow range of A4 C { where fragmentation 
occurs (for different reasons), although the timescale and collapse factor 
associated with the same M c [ is very different for different 7. 



Incorporating this approximate form for time dependence into 
our Monte Carlo approach is straightforward, since we can simply 
evolve the contribution from each mode, i.e. each A<5 ; in Eq. [5] 
according to Eq.[52] For a "timestep" At, this is equivalent to taking 



A6j(t + At) = A5j(t) exp(-Af/r) 



(54) 



+ K \/ AS ( 1 - exp (-2 At /t) ) 

w A<5j0) (1 - At /t) +11 ^2 AS At /t 

where 1Z is a Gaussian random number with unity variance, and 
the latter equality is the series expansion for small differential 
timesteps At. 

If we assume that the statistics of density fluctuations follow 
those of compressive (longitudinal) velocity differences, Eqs. [52] - 
[54] yield Lagrangian velocity structure functions (moments of the 
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Figure 12. Illustration of a "fragmentation tree," constructed as described in 
§ |1 1| for a collapsing cloud (as § |9.2) which is initially a last-crossing with 
■Mcl.o = 1 in our standard model. The cloud begins with some initial mass 
M c i at ; = 0, with no self-gravitating (fragmenting) sub-regions. Time since 
t = (in units of to from Eq. |63| approximately the cloud crossing time) is 
on the vertical axis. Each new circle ("branch" end) represents the formation 
of a new self-gravitating sub-region (last-crossing) within the cloud that 
then internally collapses and fragments, forming at the time where the circle 
is placed (and connected via the branch to its "parent" and descendant sub- 
clouds). Circle size scales with the log of each sub-region/fragment mass, 
from 0.01 M c \ (the minimum value plotted) to M c \. 



distribution of velocity differences measured for Lagrangian fluid 
parcels as a function of time sep aration) S^(A t ) = ( \\(t) — \(t + 



Af)|") oc At" /2 , identical to the |Kolmogorov| ( |1941| model (see 
Yakhot 2008). For |Kolmogorov| 1941 1 > this follows from the as- 
sumption that the statistics are strictly scale-free in the inertial 
range (independent of integral scale and viscosity); but this is 
equivalent to our derivation assuming un-correlated modesP^How- 
ever, as in the time-independent case, experiments indicate that in- 
termittency violates these assumptions and produces distinct struc- 
ture functions, which we discuss below. 

9.1.1 Generalization for Polytropic Equations of State 

It is straightforward to generalize this for a polytropic equation of 
state, using the approximation from § |3.3| that the distribution is 
locally Gaussian with AS(R, R - AR) -t AS(R, R - AR, p). We 
still need to solve the nonlinear Langevin equation, but can follow 
the same procedure as before, simply replacing steps in spatial scale 
with steps in time. 

9.1.2 Generalization for Intermittent Turbulence 

Generalizing this for intermittent turbulence as we approximate it 
in § [7] is more complicated, because the statistics follow a discrete 
Poisson distribution rather than a continuous distribution. 

Consider for simplicity the purely quantized log-Poisson case. 
Recall, when stepping from one scale to another, we have pr 2 — 
j3"p (R1/R2) 1 where m a Poisson random variable with mean 



A, and 7' is related directly to the variance (A). We can write the 
step as 



As = \n(pR 2 /p Rl ) = Aas(1 —ftp) +m\ap p 



(55) 



withP(m) = (A AjS /m!) exp(-A A s) and A as = AS/(ln/3 p ) 2 . This 
suggests that we can think of the change in density across a given 
scale for a trajectory as the cumulative result of a Poisson rate pro- 
cess. In each timestep Ar, we can generate a Poisson Am with 
mean and variance A(Af, AS) = (AS/(ln/3 p ) 2 ) (2 Af/r), and take 
m to be the sum of Am (i.e. result of the time-integral) over the cor- 
relation time, i.e. from t = t — r/2 to t (or N = r/(2At) steps)]^] 

Asj(t + At)= [— AAs(l-ft) + Am,ln/3 p ] 

1 / 2 Ar \ Am i / 2 Af \ 
P(Am,)= A^! (~ Aas ) ex P(-— A ^'J (56) 



In this particular formulation, then, the variation AS is 
contributed by some "number of events" in a Lagrangian vol- 
ume within a correlation time, with expectation value Aas = 
AS/(ln/? p ) 2 . For larger |ln/3 p | (i.e. more highly intermittent tur- 
bulence), the "number of events" is smaller, hence larger non- 
Gaussian corrections appear; while as fi p — > 1 the variance is as- 
sumed to be distributed over an infinite number of small events, 
leading (via the central limit theorem) to Gaussian statistics) 22 "] 

This form for time evolution is motivated by cascade mod- 
els (§ [7| for the intermittent Lagrangian moments and structure 
functions of the compressive velocities (S;(Af) = (\\{t) — \(t + 
A?) |")), which appear to also follow a log-Poisson hierarchy (even 
in highly compressible turbulence; see e.g. |Konstandin et aT 2012a|> 
as do the Eulerian structure functions (Dubru fle||1994| She &| 
Waymire 1995). Similar velocity scalings can be derived from sym- 
metry considerations based on the Navier-Stokes equations ( Yakhot 
|& Wanderer|2004|[Ya"khot|2008^ . Applied to the velocity moments, 
the above scaling predicts S„ comparable to numerical simula- 
tions ijKonstandin et al. 2012al. Moreover, direct observations of 
density fluctuations in the solar wind exhibit Lagrangian structure 
functions with consistent log-Poisson scalings (critically, with con- 
sistent parameters in density and velocity fluctuations, although 
in somewhat better agreement with the alternative intermittency 
model in Appendix|B] see Marsch & Tu 1994; Sorriso ^Valvo et aT] 
|1999| and references therein). This also appears broadly consis- 
tent with the evolution of scalar concentration distributions, which 
exhibit non-Gaussian features in simulations ( |Pan & Scannapieco] 
2010 1, though further study is needed(f] 



We include the factor of 2 here so that r has an identical meaning as in 
the Gaussian case: starting from As = 0, the variance for an ensemble of 
trajectories will grow in the linear regime (t — tg <C t) at the same average 
rate as in either case, and the correlation between fluctuations will decay 
with the same average correlation time. 

22 If we allow for continuous exponential damping of fluctuations (on a 
correlation time) and integrate over sufficiently large times, rather than sim- 
ply counting fluctuations in a fixed time window, we recover the thermody- 
namic intermittency models discussed in Appendix [b] which are qualita- 
tively similar but conv erge to the lognormal statistics more rapidly. 

23 As noted in § 1 8 -5 1 & [X] since our predictions depend primarily on the 
second-order and lower moments of the density field, the exact form of 
the intermittency model we adopt makes little difference provided that the 
variance (second-order Lagrangian structure function) obeys a similar linear 
scaling in time (Sj ~ Af) to our assumption. This is equivalent to extended 
self-similarity, so is true in nearly all multi-fractal models. 
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9.2 Spherical Collapse of a Barotropic Cloud 

We now consider the evolution of fields in time in a collapsing (or 
expanding) background. The specific case of interest for many ap- 
plications in gravo-turbulent fragmentation is a contracting, self- 
gravitating cloud or "clump." 

Consider a spherical, marginally self-gravitating overdensity 
or cloud, with mass Ma, radius Ra, cloud-scale averaged rms tur- 
bulent velocity V(, c i(i?ci), average sound speed at the mean cloud 
density c c i(p c i), and virial parameter Q': 



Cd(Ai) + v,.d(^ci) 2 = Cd (1 +Ma) = Q 



(57) 



where Q' ~ 1 is constant and Ma = v<, c i(^ci)/ c cl(Pci) is the cloud- 
averaged Mach number. The total cloud energy is then -GM^/Ra- 
Assume that the gas can be approximated as a polytrope c 2 s oc p 7_1 
over some limited dynamic range in p (it is trivial to allow 7 to 
change as we follow collapse). This gives 

'R A \ -3(7-1) 
- ^0 ■ 



2 
Ql 



: C 



2 --3( 7 -l) 

- Cor 



(58) 



where Ro and cq are an initial radius and (mean) sound speed of the 
cloud and r = Ra/Ro- The virial relation above gives 



l+M 2 d = Q 



, GMg 



(l+Mio)? 



1+3(7-1) 



(59) 



where M c i,o = V(, c i(i?o)/co is the initial cloud-scale Mach number. 

Let us assume that there is no source to "pump" turbulence 
other than the gravitational contraction of the cloud itself; then it is 
well-established that the kinetic energy of the turbulence decays in 
about a crossing time (dissipating in shocks or via the cascade). So 



AE, 
~dJ 



2 R c i/v,,a 

n Q n / 2 g 3/: 



(60) 



.5/2 



—5/2 
r 1 



where r\ ~ 1 is a constant calibrated from numerical simulations 
(e.g. |Kitsionas et al.|2009"l|Pan & Scarmapieco|2010| >. 

Since (by assumption) the cloud is dissipating and not gaining 
energy from outside, it must contract. If this occurs with approxi- 
mately constant virial parameter Q' , then the change in energy from 
contraction is given by 



~d7 



GMl dR A 
Rl dt 



GM L A 1 df 



dt 



Ro 



(61) 



Without external energy input, these must be equal, so (af- 
ter some simplifying algebra) we obtain the following differential 
equation for the evolution of the cloud size: 



dr 
dr 



-1-3(7-1) ,3/2 



1+^1,0 

where we define T = t/to and to is the characteristic timescale 

2 /GM c i\-i/2 



to 



Rl ) 



(62) 



(63) 



VQ' 3/2 V 

which is proportional to the initial dynamical time of the cloud. 
9.2.1 Some General Consequences 

Some general behaviors can be discerned from Eq.[62] Since f and 
Ma,o are positive and r(r = 0) = 1, dr/dt < and the cloud 
shrinks. However for 7 > 4/3, 1 — 3 (7 — 1) < 0, so the second 



term will vanish (or, if there is some overshoot, become negative 
and drive expansion) and the contraction will cease at 



^stable — 



Ro 



= (1+^5,,)^; ( 7 >-) 



(64) 



Physically, c s grows faster than v,. c i in this case, stabilizing against 
any further contraction (until this energy can dissipate) and frag- 
mentation at this 7. 

For 7 < 4/3, however, the cloud collapses to r — > in a finite 
time t ~ [A^ci,o/(l +Ml L o)~\ l/1 to. As it does so, Vt,a grows faster 
than c c i so M c i increases as r decreases. Thus either with large 
initial A4a.o 3> 1 or as f — > 0, we reach .A/fd 3> 1. In this limit 
the r 1_3( -' ) ' _1 - ) term vanishes and df/dr w — f~ 1 ^ 2 , i.e. collapse be- 
comes scale-free. In this limit, since the cloud dynamical time at 
any instant is ~ (GM a /R 3 cl )~ 1/2 oc r i/2 , and |dr/dlnrj w f i/2 , the 
cloud spends an equal number of local dynamical times at each log- 
arithmic interval in size as it collapses (even though the collapse to 
r — proceeds in finite absolute time). And in this limit, because 
the turbulent term dominates, collapse becomes independent of 7. 

For isothermal gas, we can analytically solve Eq.|62| 



r(r) 7=I =2[ x /z(l-z)- 1 -sin-'(^)] 



1 

IM' 



Vl+M> 2 ' i l-M ci 



(I+^CLO)" 



(65) 



We can trivially follow collapse for a cloud with one 7 up to 
some p c i and/or ^ c i, then switch to a different 7 and continue to 
follow it (see|Jappsen et al.||2005] |HennebelIe & Chabrier||2009[ 
|Veltchev et al.||2011| >. Thus the derivation above can be used to 
treat collapse with a complicated barotropic or bivariate equation 
of state c' s (p, R). 

9.2.2 Modification to Density Field Evolution 

As a parent cloud globally contracts, then, all scales R shrink, and 
densities correspondingly scale up. However, recall that at fixed 
A^ci the dimensionless quantities that determine the last-crossing 
distribution depend only on the relative values R/R c \ and p/pa, 
which are independent of the cloud shrinking. The only effect, then, 
on the fragmentation process in these dimensionless units in a col- 
lapsing or expanding system is through the evolution in Ma- 

So we should consider all quantities in terms of the scaled 
variables R/R c \ rather than the absolute R: essentially following 
"Lagrangian" modes k/ka =constant. Up to variation in Ma, then, 
these Lagrangian modes follow the evolution derived above for a 
stationary system, and we can trivially construct the trajectories 
S(x = R/Rd) at each time and compare them to the dimensionless 
B{x = R/Rd) = constant. 

However as the cloud contracts Ma does evolve; we therefore 
must take B = B(x, Ma) and AS = AS(x, Ma) at each time. This 
leads to time-dependence in the barrier B and variance 5 at fixed 
x — R/Ra'- from Eqs. 
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it follows that 



dBi 
drl 

d5| 

drL, 



dM 



(l+A4 d )[(2-7)/5, 



CI"! I 



p<-> 



+ b 2 M 2 A 



XP~ 



+M 2 A xP-n 
1 



dr 2 dr I 



A Mi 



(p-l)M 2 At 



(66) 



where p = p/pA and p a n = p a \t{x, Ma)/ Pa- At fixed x, this means 
it becomes easier to cross the barrier at all scales as the cloud con- 
tracts when 7 < 4/3 (and the reverse for 7 > 4/3). If Ma 3> 1, 
AB/At — s> 0, i.e. the barrier becomes self-similar in the turbulence- 
dominated regime, while the variance S systematically increases 
with A^ci as Ra decreases. 
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10 RESULTS FOR TIME-DEPENDENT FIELDS 

We now consider some basic consequences for the time-dependent 
fragmentation process in a time-evolving turbulent field. 

10.1 Global Mass Function Evolution: The Rate of 

Formation of Bound Structures and "Fragments" 

First, we consider the global evolution of the first-crossing and last- 
crossing distributions. However, care is needed. For the stochastic 
differential equations above, the instantaneous rate of change of a 
quantity like the mass function is always formally divergent (be- 
cause the discrete "events" operate as delta functions). Time evo- 
lution can only be defined over some finite interval. But averaging 
over any finite interval requires some decision regarding both the 
initial conditions of the turbulence, and what happens to regions 
that do become self-gravitating. These choices are not unique and 
will lead to different conclusions. 

10.1.1 Instantaneous "Formation Rate" from Smooth ICs 

One example is shown in Fig. [7] Here we consider the rate of 
initially forming first/last crossings of mass M, in a field that 
develops from initially uniform density. We evolve our standard 
model assuming all trajectories begin at p — po, and immediately 
count the resulting bound objects the moment the trajectory has a 
crossing. We stress that if we began with more developed turbu- 
lence/inhomogeneity, we would obtain different results, and if we 
allowed the bound objects themselves to evolve, they could become 
more massive in time (we do not allow this, simply recording their 
mass when they first become self-gravitating). 

The interesting question is how the resulting rates-of- 
formation (dn/dt) compare to the instantaneous MFs (n(M) — 
drc/dlogM) in fully-developed turbulence. For first-crossings, 
dnf/dt is flattened relative to n/(M). If we simply multiply n.f(M) 
by a constant (say, the inverse crossing time at the most unstable 
scale -Kt^mm] or, closely related, the scale where M dnf/dt is max- 
imized Rfo£), we predict a dn//df oc M~ a with slope a < 2, in 
poor agreement with the full calculation. This is not surprising: we 
should expect objects to "develop" on something like the crossing 
time t(R) on each scale, since this is the timescale for the turbu- 
lence to develop (and, formally, the correlation time over which the 
density field on that scale "resets"). If instead we compare dn//df 
to nf(M)/r(R[M]), the shapes are in much better agreement. But 
while nf(M)/r(R[M]) has the right shape, it is biased towards 
higher-mass objects than d«//df by a factor « 2. We obtain quite 
a reasonable fit (to ~ 0. 1 dex) if we shift the MF by this factor, 
i.e. dn f (M)/dt m «/(2M)/V(S[2M]). We discuss this high-mass 
truncation below. 

Based on this, we might also expect the last-crossing rate- 
of-formation to follow the "static" ni(M) divided by the local 
crossing time t(R[M]); but it does not. Instead, the rate of for- 
mation dni (M)/dt is much closer to the instantaneous MF ni(M) 
multiplied by a constant. This constant is approximately the in- 
verse crossing time at the maximum of the first-crossing rate-of- 
formation function (i.e. 1/t[R™™]). 

Why does this occur? Recall, as discussed above regarding 
the predicted correlation functions, isolated last-crossings (i.e. last- 
crossings arising independently on small scales) are very rare, since 
most of the power in density fluctuations is on large scales. This 
means that last-crossings are "seeded" by first crossings on a larger 
scale. This will occur whenever dB/dS 1. Consider the follow- 
ing (extreme) illustrative example of this. On some large scale 
Ro, the variance in density fluctuations contributed by modes with 



k ~ I/Rq is a large So 2> 1. Between this and some much smaller 
scale Ri <C Ro, the variance contributed from modes on each scale 
drops very sharply, so the total contributed is some very small 
Si <C 1 <C So. The value of the log density S — ln(p/po) at Ri is 
therefore a gaussian random variable with variance S — So + Si , the 
convolution of the contribution from large scales Ro and smaller 
scales Ri < R < Ro. But since So 2> Si, 5(Ri) is hardly altered by 
the modes on scales R<Ro- i.e. it is almost entirely set by the fluc- 
tuations around Ro. But assume that the barrier B(R) rises smoothly 
and continuously at R < Ro, so it continues to increase while S in- 
creases very weakly on smaller scales (dB/dS 3> 1). If B(R) is a 
monotonic function of R, then the "location" R where this trajec- 
tory 5 will have its last crossing just depends on 5; but the value of 
S is almost entirely determined by the fluctuations at ~ Ro - In other 
words, the number density and scale of last-crossings is set not by 
fluctuations near that scale (which evolve on the crossing time at 
that scale), but by fluctuations on much larger scales, and so are 
"seeded" on these longer timescales. 

10.1.2 The Integrated "Collapse Rate" 

An important integral quantity is the rate, integrated over the mass 
function, at which mass is collapsed into bound objects. This is 
straightforward to calculate for a given assumed initial condition 
(e.g. the constant-density ICs above) by simply integrating over the 
rate-of-formation in the MF. 

We can derive a reasonable approximation to the full nu- 
merical solution with the following simple arguments. Recall, 
the first-crossing mass is concentrated near the most unstable 
scale /ffi/min]. At a given instant near this scale, we can approx- 
imate the MF with the solution for a linear barrier (Eq. [TTJ; 
we can also note from Fig. [2] that around this scale p alt is ap- 
proximately constant (dfi/dS is small, or in the linear barrier 
So > pS). Integrating the mass density down to some 5 m ax with 
these approximations gives a total collapsed mass density F = 
Pcoiiapscd/po = erfc[So/V25n,ax] (j ust the integral over the "tail" 
of the log-normal). Using R[vmm\ ~ h (Eq. |44| l in the expres- 
sion for p cnt (Eq. 14 1 gives B Q « \n[Qk~ l {\ + k 2 )]. In Fig. [2] 
we see that the 5 max of interest (the largest S before v begins 
to rise rapidly) corresponds to the maximum in dS/dlnS (be- 
low this scale, 5 grows very slowly, so dB/dS becomes large 



and MF is suppressed). For S(R) of the form in Eq. 12 where 
we can write dS/dlnfl = ln[l + b 2 v 2 /(c 2 + k 2 Ic~ 2 )] = ln[l + 
b 2 Ml\kh\ l -"/(l+K 2 (l + M 2 h )\kh\- 2 )], this is approximately 
5 max « ln[l +0.5h 2 Mj,/(k 2 (l+Ml)) ( "- ,)/2 ]. 

Finally, since the density field is randomized on a cor- 
relation time about equal to the turbulent crossing time, 
the rate of "filling" the high-density tail of the PDF is 
~ F y<cros 8 (^[>min]) « F / t aoss (h) , where t CI0S s(h) = h/v,(h) = 
cT g [h]Q- 1 /v,[h] = {l+M 2 h ) l/2 /(M h n). This gives the approx- 
imate integrated collapse rate 

1 dM co llapsed 



M tot Sl 



dt 



(67) 



: Erfc 



a\n[Q{l+R 2 )/ii] 



s/l+M 2 L ^2 ln(l + 0.5b 2 M 2 /[k 2 (1 + M^)]( P -i)/2) J 

where the (1 + k 2 )/k term should be replaced with "1" if a "hard" 
upper barrier/injection scale is used instead of an angular momen- 
tum barrier (i.e. if k = 0), and a w 1 .0 — 2.0 is a fudge factor rep- 
resenting the detailed integration over the shape of v(R). 

In super-sonic turbulence (Mi, 3> 1), the collapse fraction is 
order-unity (tens of percent) per dynamical time, with weak depen- 
dence on Q or k. And it is a logarithmically increasing function of 
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Mi, at fixed Q, while it decreases with Mi, if Q oc Mi, and Q > 1 
(as in a <2 > 1 disk primarily supported by supersonic turbulence). 
This is discussed in more detail in Paper I; there this is shown to 
agree well with the results of a number of idealized "turbulent box" 
simulations (Fig. 1 1 therein)^ 

In sub-sonic turbulence, however, the collapse fraction is cut 
off rapidly (at fixed Q), in line with the exponentially suppressed 
MF we saw previously. Expanding the above gives a mass frac- 
tion collapsed per dynamical time of ~ Mi, exp (— v 2 /2) / (p V2ty) 
where v = M^ 1 In [2(1 + k 2 )/k]. Note, however, that if there 
were no angular momentum support ((1 + k 2 )/£ — > 1) and the disk 
were marginally stable or unstable (Q < 1), then the suppression is 
far weaker, with collapsed fraction ~ Mi, per dynamical time. The 
exponential suppression in sub-sonic turbulence relies critically on 
a thermal or angular momentum barrier sufficient to stabilize the 
disk in the "traditional" Toomre sense. 

10.2 "Top-Down" vs. "Bottom-Up" Structure Formation 

Another key general question in the time histories of forming ob- 
jects is whether bound objects evolve in "top-down" or "bottom- 
up" fashion. In "top-down" formation, objects initially form at large 
masses and then fragment into smaller objects - so much so that the 
"first-crossing" mass is typically depleted in time. In "bottom-up" 
formation, small objects assemble hierarchically into larger objects 
(this is the standard "hierarchical" structure formation in cosmol- 
ogy)- 

We examine this in Fig. [8] In order for bound objects to 
"evolve" in any sense, of course, we must allow them to have 
a finite lifetime (unlike the case in Fig. [JJ where bound objects 
are immediately removed). We therefore follow the same proce- 
dure as Fig. [7] assuming all trajectories grow from uniform initial 
p — pa, but in this case when a structure becomes bound, we allow 
it to "survive," and continue to evolve (freely allowing the density 
modes to evolve in time) for some timescale /bound equal to a fixed 
multiple of the crossing time at the first-crossing scale (fcross 
When this time expires, the object is removed and the trajectory is 
"replaced" with one at the mean density. Since the lifetimes are 
finite, we plot the time-averaged MF that results in each case. 

Unsurprisingly, when fbound -C /cross, the MF normalization is 
suppressed (since objects form on finite times but have short life- 
times), and the shape is close to the "rate of formation" in Fig. [JJ 
because objects have little time to evolve after they initially cross 
the barrier. As we increase fbound, objects are not only more common 
(being longer-lived), but allowing the turbulence to evolve within 
such objects for a longer time, we see the first-crossing MF ex- 
tend to higher "maximum" masses, while the last-crossing MF ex- 
tends to lower "minimum" masses. This continues, eventually sat- 
urating when fbound S> fcross, in which limit we recover the analytic 
prediction for the MF in fully-developed (non-time dependent) tur- 
bulence. It appears, therefore, that with increasing time for bound 
objects to evolve, structures "develop" at high masses in "bottom- 
up" fashion, and at low masses in "top-down" fashion. 

Recall for lognormal density distribution, the integrated mass 



24 Specifically, Eq. [67] appears to agree well with the results of idealized 
forced turbulent box simulations from Jv[ ~ 1 — 100 with pure hydrody- 
namics plus gravity ( Vazquez-Semadeni et al. 2003) as well as hydro, grav- 
ity, and MHD jPadoan & Nordlund|201 1^ 7 as well as galactic simulations 
with cooling, star formation, and feedback from stellar evolution I Hopkins 
|et al.|20 1 la b 2012a bl, with the simulations scattering by about a factor 
~ 2 about the predicted relation. The scatter appears to be largely related to 
the non-constant rates of collapse in the simulations. 



fraction which is self-gravitating on any given scale (independent 
of first/last-crossings) is Erfc(i/„,/\/2), where we use v m = B„,/S 
to denote the variable v with respect to the mass-weighted den- 
sity PDF (so B,„ = ln(p tT i t /po) — 5/2, as opposed to the volume- 
weighted PDF where B v = In (p cr i t /po) + 5/2). If all trajectories 
begin at u(R) ~ 0, then the dispersion in v grows at each radius 
with the crossing time t(R). So if we truncate at some fixed frac- 
tion of the mass distribution, the "upper envelope" f U p P ci grows with 
dfupper/df ~ t(R)~ . In some differential timestep df, then, this 
moves the zipper (B) crossed by the upper envelope by some Av 
corresponding to a shift AR or AM in the size/mass scale which is 
now collapsing at the threshold rate. The "collapse scale" therefore 
migrates at a rate 



dfic, 



dz/ up 



dr 



df 



dv„, \ -' 
dR 



Vt (RaM) ( dfm \ 



dR Jr, 



(68) 



If dR co \\/dt < 0, then collapse cascades "top-down": large scales 
become self-gravitating first and, with time, smaller and smaller 
scales subsequently develop self-gravitating structure within these 
parent scales. This is the typical Jeans-style fragmentation cascade. 

But if d/? C oii/df > 0, then collapse flows "bottom-up," with 
larger and larger "parent structures" collapsing subsequently. This 
corresponds, physically, to two processes: (1) accretion of larger- 
scale material by the central, already-bound clump when compres- 
sions make the "parent" region sufficiently dense, and (2) clump- 
clump mergers. Indeed, in simulations of galaxy structure that re- 
solve only the largest molecular clouds and in which those clouds 
are predominantly supported by rotation, they grow in time via 
mergers and accretion as predicted here (see e.g. |Dobbs|[2008[ 
|Tasker & Tan|2009[|Hopkins et al.|20"T2b) . 

Since df upP er/df ~ r~ is always positive, the sign of d/fqoii / df 
is the same as that of dv„,/dR. This gives a time-independent crite- 
ria for "bottom-up" growth, as a function of B and 5 alone: 



dtfcoii ^ . . ff dS IS 
—r— > iff - <-- 
df dS 2 5 



(69) 



For a linear barrier (B = Bq + pS), this is just 5 < Bo/ p. So on 
sufficiently large scales (small 5), a system with an approximately 
linear barrier grows "bottom-up," while on smaller scales, struc- 
tures grows "top-down." 

This corresponds well with what we see in Fig. [8] Below ~ 
h, S increases with decreasing R (because it measures the sum of 
contributions from fluctuations on different scales), while B also 
increases (in either the turbulence-dominated or thermal-dominated 
regimes), hence dS/dS is large and positive (see Fig.[2|. Since (as 
we show earlier) the linear-barrier approximation is not bad, it must 
be the case as 5 becomes larger that dB/dS 3> 6/(25), so structure 
grows "top-down." 

Above ~ h, we see in Fig.|5]a sharp transition where df?/d5 <C 
0. The presence of an angular momentum barrier (k) makes col- 
lapse more difficult on large scales (B(R) increases with R, while 5 
must decrease to ensure mass conservation). In this limit, dB/d5 < 
B/(25) so structure growth must proceed "bottom-up.'|^]since this 
follows the sign of dv / dR, the scale where growth transitions from 
"top-down" to "bottom-up" is just the maximal instability (quasi- 
Toomre) scale defined above. If the system is turbulent and has no 



25 In the case of cosmological structure, B constant, so the condition for 
"bottom-up" growth is always satisfied and hence traditional "hierarchical" 
structure formation emerges, even though the form of the time evolution of 
perturbations (linear growth vs. the stochastic growth here) is quite differ- 
ent. 
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angular momentum barrier, then all scales collapse "top-down" - 
this is the standard Jeans collapse-style form of fragmentation. 

10.3 Time-Dependent Fragmentation of Collapsing 
"Last-Crossing" Clouds 

Now we consider the behavior of a collapsing polytropic cloud. 
If the cloud already has "last-crossings" contained within it, each 
such crossing collapses on its own appropriate timescale, while the 
"parent" cloud continues to collapse itself. We can track each such 
collapsing object separately, so the key problem that needs to be 
addressed is the evolution of a cloud which is collapsing on its 
last-crossing scale, i.e. we can without any loss of generality focus 
on the problem of a single cloud which is a "last-crossing" within 
some parent cloud. Shortly below, we will show how multiple such 
histories can be "stitched together" into a full history for a parent 
cloud with multiple sub-clouds. 

10.3.1 Methodology 

To do so, we consider the collapsing cloud background as derived 
in § |2.4| i.e. a collapsing spherical polytropic cloud. For simplicity, 
we neglect the angular momentum (k) term, since this is usually 
negligible at the scale of last-crossings in any case. 

First, we construct a cloud with appropriate initial conditions. 
Since the largest collapsing scale without fragmenting sub-scales 
in the initial conditions is the last-crossing scale, we have by defi- 
nition p = p cr it at R = i? c i = i?i ast , i.e. Q' — 1 in § 2.4 We need to 
assign some initial conditions for the density field on scales R<R c i- 
One option would be to assume the initial density is uniform on all 
scales. But this is not consistent, strictly speaking, with clouds that 
form out of a "parent" density field. To match that case, we can sim- 
ply generate an ensemble of trajectories in the usual Monte Carlo 
fashion, chose those which have a last crossing on some scale, and 
(for now) simply discard all information on scales above the last- 
crossing scale (retaining the trajectory "below" the last crossing as 
the "initial" density field 8(R)). We can also generate an ensemble 
of these trajectories to fully sample the density PDF, by repeatedly 
drawing trajectories with the same last-crossing scale. In practice, 
it is numerically less expensive (but mathematically identical) to 
simply generate trajectories starting at a parent scale of R c \, and it- 
eratively discard and re-generate those which include subsequent 
crossings on smaller scales. 

This defines the initial density field. For a given M c \.o of the 
cloud (the mach number at the initial last-crossing/cloud scale), 
then, the cloud evolves according the equations above (§ |9.2| l. The 
absolute size shrinks, and the density field undergoes the random 
walk described above, but with Ma increasing and this introduc- 
ing corresponding evolution in the variance and barrier with time. 

As the density field evolves while the cloud collapses, even- 
tually, new crossings will appear on scales R < R i. Each of these 
should itself collapse on its own appropriate timescale; this is just a 
new collapsing "last-crossing" cloud with its own appropriate ini- 
tial Mach number. We therefore record the point at which each such 
crossing appears (along with the mass, size scale, and Mach num- 
ber of the crossing), and remove that trajectory. The surviving tra- 
jectories continue to evolve as before. 

Because, for the simple model considered here (for 7 < 4/3), 
the cloud will collapse to arbitrarily small size (A4a — > 00), and 
even though collapse proceeds in finite absolute time it spends 
an equal number of dynamical times in every logarithmic interval 
in radius, all trajectories will eventually develop a new crossing. 
We therefore evolve each last-crossing until all trajectories have 
crossed. 



Finally, note that in this scenario, for a given (fixed) 7, b, 
and /3, the only parameter that enters into the equations governing 
the cloud evolution is the initial Mach number Mci.o of the cloud 
(which, for a fixed parent disk in which the cloud forms, is equiva- 
lent to the initial size and/or mass scale of the cloud). So the "frag- 
mentation history" we calculate in this manner is a one-parameter 
family, in Mcio(Rci[t = 0]). 

10.3.2 Results 



Figs. |9|10| show the resulting fragmentation histories for last- 
crossing clouds with different initial Mach numbers, in our stan- 
dard model. 

First, we consider last-crossing clouds with initial Mach num- 
bers M c i,o = 0.1, 10, and follow their fragmentation in time. We 
specifically plot the MF of fragments produced following the cloud 
collapse, at several times. Most fragmentations initially occur near 
the "parent" cloud scale - this is the maximal instability scale 
within the collapsing cloud (in Eq.|44] the maximal instability scale 
for p = 2 and k = is « 0.7 R c \ for A4a ~ 0.1 — 1). This is con- 
sistent with both the clustering and time evolution seen thus far - 
isolated last-crossings on small scales do not typically develop "in- 
dependently," but are "seeded" by larger scales. Thus even allowing 
last-crossings to collapse evolve, and fragment in time, their char- 
acteristic last-crossing scale remains similar at each stage of col- 
lapse ("evolving" only logarithmically as each sub-cloud continues 
in its own collapse). Last-crossings on a wide range of scales, ev- 
ident in the global last-crossing MF, are a consequence of modes 
"frozen in" by fluctuations on much larger scales, rather than the 
evolution of a cloud or clouds with fixed large-scale modes (the 
isolated collapsing cloud considered here). 

Note also that when the initial Mach number is small A4 c i,o i$ 
1, it takes considerable time for fragmentation to develop. This does 
not occur until the cloud has collapsed so far that A4 C \ > 1, i.e. until 
the turbulence is super-sonic; at this point new fragmentation events 
begin to develop. For initially low Mach numbers, this means frag- 
mentation does not develop until R c i(t) <C Ra(t = 0). 

Since most subsequent crossings develop in a narrow range 
of scales, we can temporarily neglect the mass dependence and 
consider the distribution in time (or, equivalently, in cloud Mach 
number or size relative to the initial value) of crossings, shown in 
Fig. [TO] Here we clearly see the effects above: for clouds which are 
initially sub-sonic or transsonic, fragmentation does not develop 
until _M c i ~ 2 — 3; once clouds are super-sonic the process is self- 
similar, with fragmentations occurring over a factor of ~ 2 in A4 C \ 
(before all trajectories have crossed/fragmented). This corresponds 
to a factor of ~ 5 — 10 collapse in R c i, in a time about equal to the 
cloud free-fall time. 

In Fig. [TTJ we examine how this fragmentation history de- 
pends on other properties of the turbulence. We show for each case 
the Mach number dependence, but vary 7, b, and p (because these 
are defined as last-crossings at the cloud scale, Q factors out com- 
pletely). In general we see that the histories are quite similar inde- 
pendent of these choices, albeit with some interesting subtle differ- 
ences. We have also performed this calculation for different values 
of (3, but we find almost no effect here for any plausible values. 

There is no exact closed-form analytic expression for the time- 
dependence of fragmentation, but we can approximate it via the fol- 
lowing. The collapsed fraction is j dS'ff(S'); approximate f/(S) 
with the solution for a linear barrier, assume that the density distri- 
bution can reach equilibrium at each logarithmic stage of collapse, 
and approximate S by Taylor expanding near the scales R ~ R c \ 
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where fragmentation events occur. This leads to: 



2R 



d/frag 

ir~" Jut " r° dtj 

So = (3 -p)- 1 In (1 +b 2 M 2 cl ] 



(70) 



For small M c \ 1, this simplifies to approximately 
d/trag/dr ~ M~\ exp[-2b~ [ M^ 1 ] < 1 (where again r = f/fo 
defined in Eq.|63|>. Fragmentation is exponentially suppressed, be- 
cause there is simply not enough variation in the density field 
on scales below the last-crossing (A4 C \ Cl, so 5 < 1 within 
the cloud); in other words, the last-crossing scale is "frozen" and 
though it collapses, it cannot develop new self-gravitating structure 
on small scales. 

For large M c \ > 1, d/f rag /dr ~ ip~ 2 exp (— where 
ip = t [1 — 3 (7 — l)]/[2Bo (3 — /?)]; the collapse rate grows ex- 
ponentially in time then declines when most trajectories have 
already crossed. Integrating, the fragmented fraction /f rag ~ 
exp (-So [1 - t /t]) where r = 2 (3 - p) [1 - 3 (7 - l)]" 1 . So the 
entire cloud fragments within a time t = fa To that is a couple free- 
fall times. Within that time, the Mach number grows by a logarith- 
mic factor A In A^d ~ 3 — p, and the cloud collapses in size by a 
factor Aln/? c i ~ to. So shallower power spectra and stiff er equa- 
tions of state produce fragmentation over a wider range in time and 
Mach numbers, by suppressing the power in density fluctuations on 
the maximal instability scale, but the dependence is only manifest 
in a weakly dependent pre-factor (because turbulence dominates in 
this regime, and so collapse is approximately scale-free). 

11 CONSTRUCTING "FRAGMENTATION TREES" 

In Paper I, we outline how time-evolution models of the sort de- 
scribed above can be used to construct "merger and fragmenta- 
tion trees" analogous to dark matter halo merger trees in the ex- 
tended Press-Schechter formalism. This is just the procedure we 
used above in § 1 10. 1 1 to follow the development of first-crossings 
with finite lifetimes. As in § |10.1| we considered in Paper I only 
the global statistics of first-crossings, and did not allow for self- 
gravitating clouds to undergo any time-dependent evolution (con- 
traction or successive fragmentations on the last-crossing scale). 
Having further developed these models to allow for the full frag- 
mentation cascade, we now outline the methodology for construc- 
tion of a complete "fragmentation tree." 

(1) Select (or construct) the initial conditions within the 
"cloud." This may vary depending on the problem of interest. But 
they can be generated as described in § 1 1 0. 3 . 1 1 for either smooth 
initial conditions or "fully developed" turbulence within the cloud. 

(2) Evolve the system forward by one time step At (evolving 
the density PDFs per § 10 1, and calculate the new MF dN/dM — 



V c i An /AM for last-crossings below the parent cloud scale at the up- 
dated time (as § 1 10.3 .2| >. This can be done either via Monte Carlo 
methods or (if the barrier is sufficiently simple) analytic solution. 
Additional physics in the model can be applied at this stage rj 

(3) Draw a population from this MF. There are several ways to 
do this, outlined in Somerville & Kolatt ( 1999). For example, one 
can bin the entire MF in narrow mass bins and draw from a Poisson 
distribution with (N(M, M + AM)) = AM AN /AM, but one must 



26 Many possibilities for this are presented in Paper I; for example one 
could make additional assumptions about time-dependent variations in the 
equation of state, damping of the turbulent velocities, continuous accretion 
or expulsion of gas from the clouds. 



be careful in doing so to construct the "draw" so that mass is con- 
served. A simpler (although slightly less accurate) approximation 
is to assume all fragmentations are binary, for sufficiently small 
timesteps At. Draw an Mi from the MF by matching a uniform 
random variable to the cumulative N(> M) distribution, and divide 
the total mass M into M\ (which is the "fragmented" sub-unit) and 
M — Mi (the remainder of the "parent"). The latter continues on the 
same "collapse trajectory" as before the fragmentation event. How- 
ever, if Mi > M/2 (assuming the draw was from a random volume 
element, as in this case), then M — Mi < M/2 and the "residual" is 
not self-gravitating independent of the unit Mi , so the cloud should 
be treated as if no fragmentation has occurred (i.e. M — M\ is still 
bound directly to Mi , so the total collapsing mass is still M). 

If desired, however, the smaller region Mi now collapsing can 
be followed, speeding up the collapse rate, but tracking the residual 
mass M — Mi still bound but now above the last-crossing scale. This 
will then represent a subsequent, time-dependent accretion onto the 
unit Mi after it collapses. Repeating this process for all "levels" in 
a collapsing first-crossing unit, every bound parcel of gas can be 
assigned a last-crossing sub-unit to which it is most tightly bound, 
and a timescale on which it will accrete onto this sub-unit. Thus 
accretion information is implicitly included in the fragmentation 
tree. 

(4) Remove each "fragmented" sub-unit. These are now inde- 
pendent collapsing sub-clouds, each of which can be treated with 
steps (l)-(4) in the same manner as the original parent cloud. The 
mass remaining in the "original" cloud can then be used to re- 
peat steps (l)-(4) as well, until some desired threshold is reached 
(e.g. mass or size falling below some threshold, or some physical 
criteria such as a density threshold where collapse is completed 
or shut down). In sending the new clouds to (1), either new ini- 
tial conditions can be generated, or the updated density field (with 
fragmenting/last-crossing sub-regions removed) can be used as the 
initial condition for the next timestep. Note, though, that the col- 
lapse rate for the "parent" follows the same "track" as for the mass 
with which the cloud originally formed (despite the removal of the 
sub-unit), since the fragmented sub-units still exist within the par- 
ent region and contribute to its total mass. 

This can be followed for an ensemble of clouds (correspond- 
ing to e.g. the "initial" last-crossing population instantaneously pre- 
dicted for fully developed turbulence within a disk and/or more 
massive initial cloud), sampling different Monte-Carlo histories at 
a given set of initial conditions as well as the dependence on those 
conditions. 

Fortunately, for a polytropic gas with constant index 7, in a 
spherical cloud collapsing as described in § |9.2| the dimensional 
parameters of the cloud (size, mass, density, etc.) factor out and 
the statistics of the fragmentation histories form a one-dimensional 
family in the initial cloud mach number A^d.o- So the tree con- 
struction can be simplified as follows: generate a set of trees by 
iterating steps (l)-(4) for an ensemble of clouds on a grid of A'ld.o 
values; within each tree, when a fragment forms, record its own 
■Md o = M. (if&agment, teagment) (and absolute dimensional proper- 
ties) but do not follow it further. Each such fragment can then be as- 
signed the history of one of the grid ensemble with the same Md.o, 
and so on until some final collapse threshold is reached for all frag- 
ments. 

Fig. [12] shows one illustrative example fragmentation tree, 
constructed as described above, for an initially Ma,o — 1 last- 
crossing cloud, with initial conditions corresponding to fully de- 
veloped turbulence in our standard model. 
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12 DISCUSSION 

We have developed a flexible, analytic framework to understand 
fragmentation and development of self-gravitating structures in tur- 
bulent media. In Paper I-Paper III, we used the fact that the density 
distribution in isothermal, turbulent gas is approximately lognor- 
mal to show how the mathematical excursion-set formalism could 
be applied to calculate certain properties of turbulent density fluc- 
tuations. In those papers, we specifically considered the relatively 
simple case of a galactic disk with supersonic, isothermal hydrody- 
namic turbulence; nevertheless, we showed that this could explain a 
remarkable range of properties observed in giant molecular clouds, 
protostellar cores, and the ISM. 

In this paper, we generalize this in a number of important as- 
pects. We develop a framework that allows for different levels of 
rotation (and disk mass profiles), complicated gas equations of state 
(not just barotropes but multivariate equations of state as well), 
magnetic fields, turbulent power spectra, turbulent driving mech- 
anisms (e.g. ratio of compressive vs. solenoidal forcing), intermit- 
tency, non-Gaussian statistics (i.e. arbitrarily non-lognormal den- 
sity PDFs), correlated density fluctuations on different scales, and 
collapsing (or expanding) backgrounds. We also generalize the the- 
ory to follow time-dependent development of fragmentation, and 
show how to construct "fragmentation trees," analogous to cosmo- 
logical "merger trees" for the growth of dark matter halos. 

This should enable the application of this theory to a wide va- 
riety of astrophysical contexts. As discussed in Paper I & Paper II, 
there are a number of applications to study the structure, formation, 
and evolution of GMCs, protostellar cores, voids, and other features 
in the large-scale ISM. The formalism developed here allows for a 
first approximation of the effects of stellar feedback changing both 
the thermodynamic properties and turbulence within, say, collaps- 
ing GMCs and/or cores. This is critical for a detailed calculation 
of the formation of stellar clusters or the development of the stellar 
initial mass function. Following time-dependent collapse and de- 
velopment of subsequent fragmentation within cores is necessary to 
study the formation of binary and multiple stellar systems. Allow- 
ing for non-isothermal gas, magnetic fields, and Keplerian rotation 
is key to generalize to the case of proto-planetary disks and ques- 
tions of the role of turbulence in planet formation (either by direct 
collapse or by acceleration of planetesimal growth). Highly inter- 
mittent fluctuations in magnetized, barotropic media are particu- 
larly interesting for studying shocks, mixing, transport, and possi- 
bly even triggered detonations in convective stars at the late stages 
of their evolution. Many of these systems will be studied in more 
detail in future work. Our intention here is not to specify to any 
single astrophysical system but to develop a general mathematical 
framework that can be applied to a wide range of problems where 
turbulence is important. 

Here, we focus on two key scales relevant for self-gravity: 
the "first-crossing" scale, i.e. the largest scale on which objects are 
self-gravitating, and the "last-crossing" scale, the smallest scale on 
which they are self-gravitating (below which they will not frag- 
ment), hence the smallest objects into which "parent" regions will 
fragment as they collapse. We show how a wide variety of prop- 
erties of first and last-crossing "objects" can be analytically es- 
timated, including: mass functions, size-mass relations, linewidth 
(dispersion)-mass relations, correlation functions/clustering am- 
plitudes, "dynamic range" of fragmentation (the dynamic range 
over which typical sub-regions will undergo subsequent frag- 
mentations as they collapse), formation rates with time, typical 
growth/fragmentation histories, collapse and fragmentation rates, 



and nature of the "hierarchy" of structure formation (i.e. whether or 
not these objects will tend to develop "top-down" via a traditional 
fragmentation cascade or "bottom-up" via mergers with other self- 
gravitating objects). 

To first approximation, many of these properties depend 
surpisingly weakly on the properties of the medium. This is be- 
cause both gravity and turbulence (over the inertial range) are fun- 
damentally scale-free processes. As a consequence, the mass func- 
tions of first and last-crossings tend towards near-universal shapes, 
close to Schechter functions (power-laws with exponential cutoffs). 
Over the power-law range the slope is close to dn/dM oc M~ 2 , cor- 
responding to equal mass per logarithmic interval in mass, exactly 
what we would expect for a completely scale-free process. First- 
crossing mass functions are slightly more shallow (slopes closer to 
oc M~ l 8 ), because for any physically reasonable turbulent power 
spectrum there is more "power" in turbulent fluctuations on large 
scales (so small-scale high-density regions are more likely to live 
inside of larger-scale "parent" regions, rather than be isolated). This 
is closely related to a near-universal behavior of the correlation 
functions: both first and last-crossings are strongly correlated on 
small scales, and the correlation function has a nearly universal 
power-law like shape (running as £(r) oc r~ l at small scales and 
oc r~ on larger scales). The normalization of f (r) is a function 
of mass, but scales in nearly-universal fashion inversely with the 
abundance of a given population (Eq. \5\). Last-crossings, in par- 
ticular, are always strongly clustered, because most of the power 
in turbulent velocity fluctuations (hence density fluctuations) is at 
large scales - so collapsing, small-scale objects are preferentially 
formed inside larger-scale density fluctuations. Power-law corre- 
lations between mass and radius (Eq. [47] and Fig. |2j correspond- 
ing to approximately constant surface density in the turbulence- 
dominated regime), and between velocity dispersion and radius, 
naturally emerge. 

There are two characteristic scales in the problem which set 
most of the important physics. The first is the "maximal instability 
scale" (Eq. [44}. This is the characteristic scale of first-crossings, 
and the scale on which the medium is most unstable to fragmen- 
tation. This has some similarities to, but is not the same as, the 
traditional Toomre length/mass scale (in fact it scales in opposite 
fashion with respect to some parameters like the Toomre Q). In 
a turbulent disk, this is near the disk scale height, or in a driven 
turbulent "box" near the driving scale. Above this scale, density 
fluctuations are suppressed by mass conservation and so the num- 
ber density of collapsing objects is exponentially suppressed. If the 
presence of this scale owes to an angular momentum barrier (as in 
a rotating disk), then above this barrier the sense of structure for- 
mation is also reversed from "standard" fragmentation cascades: 
more massive objects are more likely to form via the hierarchical 
merger of incompletely collapsed (angular momentum-supported) 
clouds, rather than form or fragment out of spontaneous density 
fluctuations or larger objects. For example, many simulations of 
GMC formation in galaxies, which can often only resolve the most 
massive objects, see mergers dominating the formation of such sys- 
tems, while simulations with smaller mass or force resolution see 
far fewer events (compare e.g. Dobbs 2008 ; Hop kins et al.|2012b"| >. 

The second key scale is the sonic scale (Eq. [45}, the spatial 
scale below which the rms compressive Mach number becomes 
sub-sonic (and corresponding mass scale defined by the minimum 
self-gravitating mass on that spatial scale). This has some dimen- 
sional terms in common with the traditional Jeans mass, but it is not 
the same. This defines the characteristic scale of the last-crossing 
distribution. Below this scale, density fluctuations are suppressed 
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by thermal and magnetic pressure, so once again the number den- 
sity of collapsing objects is exponentially suppressed. We stress 
that this suppression is not because thermal and magnetic pressure 
suddenly become much "better" at resisting gravity on these scales. 
Rather, it is the suppression of density fluctuations, with sub-sonic 
Mach numbers, that is critical. 

Together, these two scales define the "dynamic range" of frag- 
mentation (Fig. [4|. While most of the volume, under any condi- 
tions, tends to be non self-gravitating, most of the mass in self- 
gravitating objects is embedded in "parent" objects on the maximal 
instability scale and undergoes a fragmentation cascade down to 
the sonic scale. 

The most important effects of varying the global turbulent 
power spectrum (spectral slope, Mach number, compressive ver- 
sus solenoidal forcing), global disk (or "parent" cloud) stability 
(Toomre Q or virial parameter), and gas equation of state, are cap- 
tured in their effects on shifting the maximal instability scale and 
sonic scale. In fact, most of the differences in the predicted mass 
functions (Fig. |3j and correlation functions (Fig. |5j with different 
model parameters can be eliminated if we simply "renormalize" the 
spatial or mass scales ("stretching" or "compressing" each in scale 
to match the maximal instability scale and sonic scale in each case). 

This gives a simple intuition for the effects of most parame- 
ter variations. Higher Mach numbers produce fragmentation over a 
wider range of scales because the turbulence is super-sonic (able 
to produce large density fluctuations) down to smaller scales; if 
the global Mach numbers are sub-sonic, the total mass fraction in- 
volved in fragmentation at all becomes exponentially suppressed. 
A shallower turbulent power spectrum, anchored to the same large- 
scale Mach number, implies more power at small scales (a smaller 
sonic scale) and so promotes fragmentation to smaller scales as 
well. Increasing the global stability of the system (i.e. raising the 
"barrier" required for a fluctuation to collapse) leads to a narrower 
range around the maximally unstable scale on which collapse is 
likely, and suppresses the overall mass fraction involved in col- 
lapse. Since it is the compressive component of the Mach number 
which drives density fluctuations, changing the balance of com- 
pressive (longitudinal) vs. solenoidal modes is nearly degenerate 
with changes to the Mach number. Changing the gas equation of 
state, while it introduces non-Gaussian density PDFs and compli- 
cated correlations between modes on different scales, has the most 
important effect of shifting the sonic scale (i.e. fragmentation is 
promoted over a wider dynamic range with a softer equation of 
state, simply because the cascade can proceed further before hit- 
ting the sonic scale). In the super-sonic scale regime, the thermal 
physics of the gas has almost no effect. 

For these reasons, phenomena such as intermittency, non- 
Gaussian statistics, correlations between turbulent fluctuations 
on different scales, and anisotropic collapse are predicted to 
have suprisingly little impact on the statistics of fragmentation 
(Figs. [Al|A4[ l. To the extent that there are some effects, they are 
largely degenerate with smaller variations in the Mach number, 
equation of state, or global stability of the system. 

We show that fragmentation develops globally (in an ini- 
tially smooth system) starting near the maximal instability scale. 
First-crossings tend to form on the turbulent crossing time at each 
scale (Fig.|7J. If objects have finite subsequent lifetimes (Fig. [8}, 
then a fragmentation cascade proceeds from the largest to the 
smallest scales, with the dynamic range of self-gravitating object 
sizes/masses increasing with that lifetime. 

Last-crossings, on the other hand, tend to be born "fully 
formed" when their "parent" region undergoes a density fluctua- 



tion that pushes it above first crossing. As noted above, because 
the variance on the characteristic scale of last-crossings (the sonic 
scale) is small (by definition), they only very rarely arise via an 
entirely local (small-scale) fluctuation, but rather "ride on" larger- 
scale, larger-amplitude density fluctuations (i.e. first crossings). As 
such, they are not formed on the local crossing or dynamical time, 
but rather "seeded" by first crossings, and so form at a rate regu- 
lated by larger-scale fluctuations in the turbulent medium (the char- 
acteristic timescale of the maximal instability scale, rather than the 
sonic scale). In other words, the mass spectrum of fragmentation in 
objects at or below the sonic scale is "frozen in" by fluctuations on 
larger scales. 

We show how to follow the development of subsequent fluc- 
tuations in these objects as they collapse. Recall, fragmentation 
which is "pre-seeded" on all scales is captured in the statistics de- 
scribed above; but it is possible that a last-crossing region (which at 
a given instant contains no fragmenting sub-regions) could develop 
such sub-regions later in time as it collapses. Formally speaking, for 
"soft" equations of state (7 < 4/3), if collapse proceeds infinitely 
(to zero size) at constant virial parameter, the system must even- 
tually fragment. However, even under these conditions, collaps- 
ing last-crossing objects which are initially at or below the sonic 
scale (i.e. have cloud-scale compressive Mach numbers < 1) de- 
velop such fragmentation very slowly (Figs. |9jTTJ. In fact, signif- 
icant sub-fragmentation only occurs when - as a consequence of 
these assumptions (in particular collapse at constant virial parame- 
ter) - the cloud Mach number becomes super-sonic (A4a ~ 2 — 3). 
If some other physics eventually enters to prevent this from oc- 
curring - for example, if the equation of state becomes stiffer as 
the system collapses, or the collapse proceeds directly rather than 
converting all energy to turbulent motions (i.e. the virial param- 
eter is not constant during collapse), then collapse will proceed to 
r — > within such an "initial last-crossing" without ever developing 
significant sub-fragmentation. In this case, the last-crossing distri- 
bution predicted at a given instant will indeed represent the final 
collapsed-fragment distribution. Thus the statistics of fragmenta- 
tion, and quantities like the last-crossing distribution, are "frozen 
in" at these scales, and only weakly modified by the collapse pro- 
cess itself. In contrast, a cloud which is initially highly supersonic 
(but somehow contrived to have no self-gravitating sub-regions) 
will develop sub-fragmentation very rapidly (in one crossing time) 
as it began to collapse, independent of the turbulent power spec- 
trum or gas equation of state. 

Our intention here is to develop and present a robust, gen- 
eral framework to understand "gravo-turbulent" fragmentation and 
structure formation. In future work, we will apply the results here to 
many specific observed systems of astrophysical interest. It is also 
important to test many of the assumptions and analytic calculations 
in the model here by comparing to fully non-linear numerical simu- 
lations and calculations, especially idealized cases (for example, an 
isothermal driven turbulent box with self-gravity), in which the rel- 
evant input parameters of the model (the turbulent power spectrum, 
gas equation of state, etc.) can be completely specified. Some par- 
ticularly important questions include how the density power spec- 
trum relates to the velocity power spectrum, how this is altered in 
the presence of rotation and magnetic fields, and how intermittency 
is manifest in the log-density field (as opposed to just the velocity 
field). Many of these are quite interesting questions in their own 
right, which should inform our general understanding of compress- 
ible turbulence. In a companion paper (in preparation), we will 
present a preliminary series of such comparisons useful both for 
developing a deeper understanding of many behaviors seen in sim- 
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ulations, as well as testing and calibrating the theoretical assump- 
tions necessary here. 
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Figure Al. Density PDFs and structure functions for the polytropic (7 ^ 1) 
and intermittent (/3 ^ 1) models in Figs. |3|A4] Top: Density PDFs for dif- 
ferent polytropic index 7 (§|3j, normalized to the same Mach number at po 
(such that the variance 5(R, po) = 1)- Non-isothermal gas leads to power 
law-like tails at low (7 > 1) or high (7 < 1) densities and steeper sup- 
pression of opposite fluctuations. Middle: Density PDFs at fixed S = 1 for 
isothermal (7 = 1) gas with no intermittency (j3 = 0p = 1) and (probably 
unphysically) strong intermittency according to the quantized log-Poisson 
|She & L eveque 1 1994} model (J3 P = 0.2). The latter is skewed to low densi- 
ties, discretized, and vanishes entirely for densities above a maximum near 
~ Po. We compare the same (5 P = 0.2 model convolved with small Gaus- 
sian scatter such that 10% of S is Gaussian. This is still highly non-Gaussian 
but removes the most severe discreteness effects. Bottom: Velocity structure 
functions S„ = ( 1 8V \ ) oc , for different choices of /3 = 0p . Kolmogorov 
turbulence with lognormal statistics is (3 = 1 . 



APPENDIX A: DETAILED EFFECTS OF 
NON-GAUSSIANITY, INTERMITTENCY, AND 
CORRELATED TURBULENT FLUCTUATIONS 

Al Overview 

In two places in the text we have considered non-Gaussian, inher- 
ently correlated structures in turbulence: first, for non-isothermal 
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Figure A2. First & last-crossing MFs in non-isothermal gas (7 = 2/3 and 
7 = 4/3). We show the exact Monte Carlo numerical solution (§ |2.1| his- 
tograms), accounting fully for a non-Gaussian density distribution and (nec- 
essarily) correlated modes in the density field between different scales. We 
compare the analytic result using Eqs. |4|7| derived (for 7 = 1) assuming 
locally Gaussian statistics and fully un-correlated mode structure between 
scales, but using the appropriate p cr i t (R) for 7 ^ 1, and replacing the vari- 
ance S(R) (for 7 = 1) with S(R, p c rit) ( see Fig[2}- The two agree well. 
Although changing 7 can significantly change the MF, the dominant ef- 
fect is not the correlation structure of the density field nor the local non- 
Gaussianity in the density statistics. 



7 / 1 equations of state (§ [3j, and second, the case with inter- 
mittency, /3 P 7^ 1 (§[7p. Both of these produce non-Gaussian den- 
sity PDFs and non-normal statistics. In fact, for certain choices the 
density PDFs resemble one another - however we stress that the 
intrinsic correlation structures of the turbulence implied are quite 
different. 

Recall, in taking a "random walk" with 7 7^ 1 to sample the 
density field, we must allow for a density-dependent local disper- 
sion S = S(R, p) and solve each trajectory as a correlated random 
walk - in other words, the "steps" depend on the absolute position 
at each step and the position at every other step. With intermittency, 
in the context of the She & Leveque ( 19941 model, it is still true (in- 
deed it is an assumption of the model) that the "stepping" is inde- 
pendent of position (in log space); however the statistics are highly 
non-Gaussian and spatially correlated by the skew in the parameter 
P P and the discreteness of the Poisson distribution. 

Some effects of different assumptions about intermittency 
are illustrated in Fig. |A1| First, we compare the density PDF, at 
(for convenience) a fixed total variance 5=1. Non-intermittent 
isothermal turbulence produces a lognormal distribution with me- 
dian —S/2. We compare this to the density PDF resulting from 
strictly quantized log-Poisson statistics - the simplest mathemat- 
ically allowed PDF that matc hes the|She & Leveque| ( |1994| struc- 
ture functions - as assumed in |She & Wa ymire ( 1995|> and|Dubrulle 
( |1994} (but applied to the density statistics as described in § |7J. In 
this case the density PDF is discretized (rather than continuous), 
with "jumps" quantized in units of ln/3 p . It is also highly skewed, 
with a large power law-like tail to low densities and a sharp, ab- 
solute cutoff at modest p > po- Clearly, the discreteness effect is 
not realistic (it is an artifact of our simplifying assumptions); recall 
from § [7] that more generally the statistics in intermittent turbu- 
lence should be log-Poisson P(m) convolved with some function 
G«(ln {p/ po) /ln/3 p — m) reflecting the variation in strength of dis- 
sipative events. We therefore also consider the same model, but as- 
suming this convolution function Gr is Gaussian, with a dispersion 
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Figure A3. First & last-crossing MFs in isothermal turbulence with dif- 
ferent levels of intermittency (§|Aj. Top: We show the exact Monte Carlo 
distribution accounting for the non-Gaussian density distribution and cor- 
related modes between different scales corresponding to our "standard 
model" but different /3 = fi p and 7' parameterizing the intermittency (ac- 
cording to the quantized log-Poisson She & Leveque 1 1994 1 model for in- 
termittency). j3 = — 1 parameterizes the strength of intermittency and 
correlation between fluctuations: /3 = 1 is non-intermittent, uncorrected 
modes with Gaussian statistics; /3 — > gives fully intermittent, perfectly 
correlated fluctuations on different scales with infinitely skew statistics. 
/J = 2/3 is the "standard" intermittency seen in experiments and sim- 
ulations of isotropic hydrodynamic turbulence. /3 = 1/3 corresponds to 
more singular turbulence with infinitely thin sheets and shocks or strongly 
magnetically-dominated collisionless media. Smaller /3 are generally not 
seen but are shown for comparison. Strong intermittency - as parameterized 
in the quantized log-Poisson model - creates a sharp (but possibly artificial) 
discrete cutoff in the density distribution above p ~ po, truncating the MF 
at low/high masses. The effect is similar to increasing Q. Bottom: Same, 
but using "inverse /3" models which represent identical levels of intermit- 
tency, but reverse the sign of the skewness; this biases the PDF towards high 
densities and the MFs to a broader range, but the effect is weaker. 



equal to 10% of the total variance (the other 90% coming from 
the log-Poisson distribution). Physically, this corresponds to 10% 
the variance coming from variation in the "level of dissipation" as- 
sociated with a given structure (e.g. the size and/or strength of a 
shock or vortex), while the remaining 90% comes from variance in 
the "number of structures" in some Lagrangian volume. This con- 
volution is sufficient to make the density PDF continuous, though 
it preserves the multiple peaks, skewness, and sharp high-density 
cutoff. 

Fig. |Al| also shows the results of intermittency on the velocity 
structure functions, which are well-studied. Recall, f3 — 1 corre- 



ct 
O 

c 

"D 



CD 
O 



0.0 
-0.5 
-1.0 



(3=1 (Gaussian) 
(3=2/3 (She & Leveque) 
(3=1/3 (Boldyrev) 
-(3=1/10 

(3=1/100 



CD 
O 

c 

"D 



CD 
O 




-6 -4 -2 
log(M) [p h 3 ] 



Figure A4. As Fig. |A3| but assuming 10% of the variance in the density is 
contributed by normal fluctuations, while the remaining 90% comes from 
quantized log-Poisson fluctuations with the labeled j3. While the qualita- 
tive effect is the same (narrower MFs with lower (3, as the density PDF is 
skewed towards progressively lower densities), the dependence on f3 is sub- 
stantially weaker. The most severe effects in Fig. |A3| arise from the strictly 
discrete density PDF in the simple log-Poisson model; these are eliminated 
by even a small amount of variation which allows both positive and negative 
fluctuations. The skewness is almost identical (within ~ 2%) in each case 
to the "pure quantized" model, so non-Gaussianity alone is not the domi- 
nant effect. Here, the effect is almost degenerate with small changes in A4/, 
(from ~ 6 - 10 for /3 < 1) or Q (from 1 - 2). 



sponds to our standard (Gaussian/non-intermittent) assumption. In 
I She & Leve que ( 19941, the authors argue that the symmetries in the 
Navier-Stokes equations lead to an analytically predicted j3 = 2/3, 
for isotropic, hydrodynamic turbulence with filamentary shocks; 
this predicts structure functions in good agreement with a wide 
range of experiments (including pipe flows with/without bound- 
aries, longitudinal and transverse shear flows, buoyancy-driven tur- 
bulence, Taylor-Couette flows, and both two and three-dimensional 
rotating turbulence) as well as numerical simulations (for a re- 
view, see e.g. She & Zhang 2009 , and references therein). Boldyrev 
(12002b used the same model to derive a scaling with /? = 1/3, for 
highly super-sonic (Burgers) turbulence, with singular sheet-like 
shocks as the primary dissipative structures; this appears to agree 
with measurements of turbulence in molecular clouds (Boldyrev 
|et al.|20 02 ) and numerical simulations of other interstellar medium 
processes ( |Pan et aL] [2009 , again for a review see |She & Zhang 
2009). This also happens to be the result for highly-magnetized 
MHD turbulence (Miiller & Biskamp 2000) and agrees with the 
structure functions observed in the solar wind and corona. For 
the sake of comparison, we consider the more extreme choices of 
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P = 1/10 and 1/100. The former may correspond to some situa- 
tions where large, non-hydrodynamic motions create density voids 
that fill the large majority of the volume (e.g. the inter-galactic 
medium; see |Liu & Fa ng 2008), though this is clearly not a "tur- 
bulent" medium in the traditional sense. The latter is purely il- 
lustrative, as (to our knowledge) no turbulent self-gravitating fluid 
systems exhibit such extreme behavior. In each case, we note that 
while formally our model with Gr containing ~ 10% of the vari- 
ance produces small corrections to the structure functions, they are 
generally negligible. 

In either case (/3 / 1 or 7 / 1), it is important to examine 
how the non-Gaussianity and different inherent correlation struc- 
tures change our results. 

A2 Effects on the Mass Functions 

A2.1 Non-Isothermal Equations of State 

We saw in § |8.2| (Fig.|3} that changing 7 7^ 1 has a significant effect 
on the MF. However, this could entirely owe to the effects on the 
barrier (changing the density threshold for collapse), or the broad- 
ening/shrinking of the density distribution (i.e. the non-Gaussian 
PDF simply "looking like" a Gaussian PDF with a different vari- 
ance 5 near the densities of interest), rather than the effects of statis- 
tics being locally non-Gaussian and/or the effect of the modes on 
different scales being inherently correlated. 

To test this, we examine in Fig. |A2| the first and last-crossing 
MFs for two choices 7 < 1 and 7 > 1 . We first show the full numer- 
ical result (from Fig.|3j, using our Monte Carlo method (§[33} that 
includes the inherent correlation structure of the modes on different 
scales, their dependence on density, and the non-Gaussianity in the 
statistics. 

We compare this to the prediction if we just apply a slightly 
modified form of our analytic solution for first and last-crossing 
(Eq.|4]&|7j, derived assuming locally Gaussian statistics and com- 
pletely un-correlated fluctuations on different scales. In the ana- 
lytic solution, we simply modify the barrier B to account for the 
appropriate critical density for collapse (p C tit(R) — > Pcrii(^| 7)) and 
replace the variance 5 appearing in the equation is with the "ef- 
fective" S(R) — > S(R, pcrit) near the critical density (as defined in 

We see that this captures essentially all of the key informa- 
tion, and reproduces the full numerical solution quite well. In other 
words, the intrinsic correlation structure of modes has a quite small 
effect on the MF. Non-Gaussianity, to the extent that it matters, is 
only important insofar as it makes the density distribution more or 
less broad at the densities of interest - re-mapping to a Gaussian 
with some different dispersion gives similar results. 

A2.2 Intermittency 

In Figs. |A3|A4| we examine the effects of intermittency (as approx- 
imated in § |7j on the predicted MF, for various values of f3 = /3 p , 
and both the "standard" [She & Waymire] ( [T995) ; |Dubmlle|ft9j?4) 
formulation (which produces a distribution skew towards excess 
low-density gas) and a mirrored "inverse beta" formulation de- 
scribed below (reversing the skew). In Fig. |A3| we consider strictly 
quantized log-Poisson statistics (so the density "jumps" are dis- 
crete, as described in § |7.1[ >. In Fig. |A4| we compare the same mod- 
els but making the density PDF continuous by assuming 10% of 

27 We do have to be careful to modify the integration to avoid where 
AS(R, p cl it) between two steps becomes negative; in which case no objects 
should appear on this scale. 



the total variance is associated with a Gaussian convolution func- 
tion Gr (as in Fig.jATJ. 

Figs. |A3|A4| show that "realistic" levels of intermittency have 
weak effects on our results. At /3 — 1/10 (/3 P = 0.46), we begin to 
see larger effects, truncating the MF at the highest-mass scales (R > 
h) and shifting the last-crossing by a factor ~ 2, and by 8 = 1/100 
ifip = 0.21) we see the MF is highly compressed. In these models, 
lowering 8 skews the density distribution towards lower values and 
so "tightens" the MF range; this is very similar to decreasing the 
Mach number of increasing Toomre Q (see Fig. [3}. 

Even at extreme 8 values, much of the effect seen in Fig. [A3] 
depends on the strict discreteness of the distribution; when we con- 
sider in Fig. |A4| the same models with just a small fraction of 
the variance in a smooth Gaussian component, the differences are 
significantly smaller. Thus, for "typical" intermittency models, the 
most important quantity determining the MFs appears to be the to- 
tal variance in the density distribution, not the detailed form of the 
multi-fractal cascade model. 

A2.3 "Inverse Beta " Models 

As shown in Fig.jATJ the model from § |7.1| for < B p < 1, pro- 
duces a density PDF skewed towards low densities (because the 
distribution of m is not symmetric). However, there are physical 
situations where the skew might be reversed) 28 ] 

It is trivial to show that the skewness of the density PDF model 
will be reversed, while the degree of non-Gaussianity in log-space 
is preserved, if we simply take f3 p — > 1//3 P (or ln/? p — > — \n/3 p ), 
i.e. fi p > 1. For fixed variance S(R), this reverses the sign of 7', 
but otherwise the equations in § |7.1| are well-behaved. We denote 
these "inverse beta" models. In this case, a larger value of f3 p (since 
this model corresponds to its "mirror" with 1 //3 P ) represents larger 
deviations from Gaussianity. We stress that we are not arguing for a 
particular physical interpretation of these models, only noting that 
they are provide a potentially useful illustration of how certain re- 
sults depend on the skew of the density PDF. 

We show the predicted MF from these models in Figs. |A3|A4| 
As expected, they produce deviations from the strictly log-normal 
case with the opposite sense of the "standard" intermittency models 
{j3 p < 1); but the magnitude of the effects are somewhat smaller. 

A3 Effects on Correlation Functions 

We might expect to see a larger effect from non-Gaussian statistics 
in the correlation functions. But in Figs. |5|6| we actually see very 
little dependence of f (r) on either the equation of state 7 7^ 1 or 
intermittency parameter j3 7^ 1. And to the extent that dependen- 
cies appear in Fig. [6]at intermediate scales and some masses, they 
appear consistent with the near-universal dependence of clustering 
amplitude on number density from Eq.|51| 

This is rather surprising at first: recall the density fluctua- 
tions on all scales are intrinsically correlated when 7 / 1, and for 
ft 7^ 1 , the structure functions (higher-order moments of the veloc- 
ity/density field) deviate substantially from that predicted by self- 
similar, Gaussian statistics. In Fig. | A 1 1 we show these, in fact, to 
highlight how large the deviation actually is. 

Why then do the correlation functions not change? Recall, as 
discussed in § |8.5| most of the correlation function structure in 

28 Although this may not be related to actual intermittency: for example 
when 7 < 1, or when self-gravitating regions are allowed to collapse but still 
included in the density PDF ( Vazquez-Semadeni & Garcia 2001 , Bournaud 
|et al.|2010||Ballesteros-Paredes et al.|2011| . 
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Figure AS. Mass functions in our standard model with different degrees 
of intermittency (as Fig. [A3), but using the alternative, continuous "ther- 
modynamic" intermittency model from Castaing 1 1 996 1 described in Ap- 
pendix|E] This provides a better fit to the density PDFs in simulations. We 
consider three values of the parameter T described in the text, T = (no 
intermittency, our standard Gaussian statistics case), T = 0.05 (which gives 
identical structure functions to the ft = 2/3 She & Leveque ( 1994) m odel), 
and T = 0. 12 (identical structure functions to the ft = 1 /3 |Boldyrev|{2002) 
model). The results are very similar to those in Fig. |A3| with slightly smaller 
deviations from the Gaussian case for the same "degree of intermittency." 



Fig. [5] is not be a consequence of inherent correlations in the tur- 
bulence, because in non-intermittent, isothermal turbulence we as- 
sume fluctuations on all scales are uncorrelated. But the correlation 
functions reflect the run of variance with scale, simply as a basic 
consequence of statistics: if the run of variance with scale is weak 
on small scales, then the probability of "isolated" strong fluctua- 
tions (large positive overdensities) arising on these small scales is 
low, and most such fluctuations will be embedded in larger-scale 
fluctuations. For a given run in S(r), the correlation function is, to 
first order, fixed. Inherent correlation structures in the turbulence 
produce only second-order corrections to this in any of the models 
considered here. 



A4 Summary 

In short, the predicted quantities in this work - the width of the 
density PDF as a function of scale, and two-point clustering - de- 
pend primarily on the lowest-order structure functions of turbu- 
lence (second-order and below). These are relatively insensitive 
to the observed correlation structure and non-Gaussian statistics 
measured in intermittency (in either experiment or simulations; see 
Fig. |Al) . At realistic levels of intermittency and non-Gaussianity, it 
therefore seems unlikely that these effects will dramatically change 
our conclusions. 

We should expect that the differences between models with 
different ft and/or other multi-scale correlated structures will be 
more pronounced if we considered higher-order correlation func- 
tions (the three-point function, etc.). As shown in Fig. |Al| the struc- 
ture functions at n = 2 (the order of the two-point correlation func- 
tion) differ relatively little even for large differences in f3, but the 
differences grow with increasing order n. It requires going to much 
higher order n>5 before large differences become apparent. 



APPENDIX B: ALTERNATIVE MODELS FOR 
INTERMITTENCY 

As discussed in § [7] a quantized intermittency model for density 
"steps," while capturing intermittency in the structure functions, 
gives an unphysically discrete density distribution. 

A different model of intermittency, with continuous variations, 
is the "thermodynamic" model proposed in |Castaing 1 1996). In 
that work, the assumptions are quite different from|She & Lev-] 
|eque| |l994). However, as shown by |He et al.| ( |1998| >, this is just 
a more general consistent formulation of the same hierarchy with 
a different choice of Gr. In this model, the predicted longitudinal 
(compressive) velocity increments scale as 



P(u)du = y~] 



A e u 



(m-1) 



- dw 



(Bl) 



= h (2 v Am) exp[— (A + u)] y — dw 

where u = T~ { j In (o>/0l)| with T an "intermittency" constant and 
a, the characteristic velocity amplitude (such that (|5v|") oc a"); 
<jl is this amplitude at the driving scale. 1\ is the modified Bessel 
function of the first kind. 

The variance in u is just 2 A, with the corresponding variance 
in | In (or/a L ) | equal to 2 XT 2 . For large- A and/or large total vari- 
ance, this approaches a Gaussian in u with mean (u) = A and vari- 
ance 2 A (i.e. a lognormal in (ov/ctl) with variance 5 and mean 
S/2). 

But it is straightforward to see that this is just the general form 
of the log-Poisson statistics with 



G R (Y\m)dY 



"du 



(m-1)! 

/m 



(B2) 



where T — > | In /3 1 w — ln(/3)/9. The predicted structure 

functions scale as 

n 1 + 3T 



(B3) 



3 1+nT 

which are nearly identical to the She & Leveque ( 1994) structure 
functions for T ^ -{j' /6) In/3, at least to n ~ 20]^] 

If we follow the assumptions in §|7] that the statistics of log- 
density fluctuations follow the statistics of longitudinal velocity 
fluctuations, this leads to the volumetric density PDF: 



PQnp) din p = h{2 V Aw) exp [-(A + w)] 

\m —A m— ' 

A e W 



dw 



(B4) 



A 



(m- 1) 



- du 



l+T 



lnp 



(w>0) 



with A = \(R) = S(R)/(2T 1 ). 

Because this distribution is infinitely divisible, incorporat- 
ing it into our Monte-Carlo trajectory calculations is a straight- 
forward extension of the model in § |3.3| & |7.1| In each "step" 



29 In|Hopkins 1 2012a) we show that T = 0.05 gives nearly identical struc- 
ture functions to the standard ■y' = 2/3, ft = 2/3 quantized log-Poisson 
model in She & Leveque ( 1994); and T = 0.12 gives structure functions 
identical to the 7' = 2/3, ft = l/3|Boldyrev| 



(20021 model. 
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R—¥ R — AR with change in variance AS, we simply draw the dif- 
ferential change in density Alnp from the distribution in Eq. |B4| 
with A = AA = AS/(2T 2 ) (instead of a Gaussian). Physically, 
this has a simple interpretation in terms of the density and veloc- 
ity fields as continuous-time multiplicative random relaxation pro- 
cesses. The variable A represents some "number of events" while 
Gr represents a convolution over a Poisson waiting time distribu- 
tion for each number of events. In other words, Eq. |B4| is a station- 
ary result for exponentially damped perturbations driven by dis- 
crete random events with a constant average rate; where T~ is a 
dimensionless "rate parameter" (higher- T corresponding to rarer, 
more highly intermittent "events"). 

In [Hopkins (2012a), we specifically show that this form for 
the density PDF provides an excellent fit to numerical simulations 
spanning M ~ 0. 1 — 15, with a variety of forcing schemes, numer- 
ical methods, and magnetic field strengths. Moreover, the same T 
fit from |B4| to the density PDF appears consistent with the values 
fit from Eq. |B3| to the longitudinal velocity fluctuations, suggest- 
ing that our simple assumption of mapping between the compres- 
sive velocity and density statistics may be plausible (at least for the 
lower-order statistics that matter here). And this form for density 
fluctuations is also favored by direct observations of the solar wind 
l |Forman & Burlaga|2003||Leubner & V6ros|2005b"la) . 

For now we simply consider whether this form for the den- 
sity PDF changes our conclusions. Fig. |A5| plots the MFs for our 
standard model, allowing for different levels of intermittency (as 
in Fig. [A3}, but with the density PDF model above. For T = 0.05 
and more highly intermittent T = 0.12, the resulting MFs show 
very weak deviations from the no-intermittency (T — > 0) case. The 
sense of the deviations is identical to the corresponding quantized 
log-Poisson model, albeit smaller. So we conclude that the detailed 
form of the multi-fractal model for intermittency makes little dif- 
ference to our conclusions here. 



to S and performing some simplifying algebra gives: 



APPENDIX C: GENERAL FIRST/LAST CROSSING 
SOLUTIONS IN CONTINUOUS TURBULENCE 

We now present the general solutions for first and last crossing for 
any infinitely divisible and continuous PDF. 

First, consider first-crossing; the derivation follows Zhang & 
|HuiH200"6) , but without simplifying by assuming Gaussian statis- 
tics. Again consider trajectories integrated from 5 = (R — > oo) and 
let 5 = In (p/po), and as before consider // and n, the probability 
of a given 5 < B(S) without having a previous crossing: 



rS rB(S) 

1= / f f (s')ds'+ / n(s\s)ds (ci) 

JO J -oo 

U(S\S) = P (S\S) - [ dS'f f {S')P w (S[S]\B'[S']) (C2) 
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g 2 (S,S') = ^P (B-B'\S-S')+ 8i (B-B'\S-S') 



gi{Ay) 



d6-P (5\y) 
, oy 



MS) 



a = l-Limitjy dSP (8-B'[S'] S-S')] ., 

For the distributions here a depends simply on the symmetry 
properties of Po; for Po symmetric in 5 (e.g. Gaussian statistics), 
a — 1 /2; for Po which is one-sided, a = (the /3 < 1 or "normal" 
intermittency models) or a = 1 (the B > 1 or "inverse" intermit- 
tency models). 

For any continuous function Po(S \ S) and B(S), the functions 
gi, gz, g3, and a are also continuous, and can be tabulated so that 
Eq. |C3| (a Volterra integral equation of the second kind) can be 
solved by standard numerical methods. 

Now, consider last-crossing. We begin the walk at some suffi- 
ciently large 5, as R — > such that the probability of crossing van- 
ishes, and proceed in "reverse" direction. The integral constraint 
is 



1 



S, rB{S) 

fe(S')dS'+ / n(<5|S)d<5 



(C4) 



U(6\S) = P Q (S\S) - J dS' f e {S')P 0l (SlS]\B'[S']) (C5) 

where Poi (S[S] | B'[S']) represents the probability of having 8 at S 
given an earlier value B'[S'], but for a step in the opposite direc- 
tion from the first-crossing case. This is related to Pio by Bayes's 
theorem: 



P 0l (5[S]\B'[S']) = Pio(B'[S']\S[S]) 



Po(S\S) 



(C6) 



= P (B' - S\S' - S) 



Po(B'\S') 
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Po(B'\S') 

Again taking the derivative of both sides and simplifying, we obtain 

ril ds'Ms')g 2 (s,s') =gi(S)-&fi(S) (C7) 



i: 



dB 



g 2 (S,S')= —Poi(B[S]\B'[S'])+ / d<5 — P 01 (6[S]\B'[S']) 



BIS] 



d 



a = l — Limit 



l(S) 



dSP 0l (5[S]\B'[S']) 



with gi the same as in the first-crossing case and Poi defined above. 
From this form, again, for any continuous Pq(8\S) and B(S), the 
input functions are straightforward to tabulate and fi (5) can be de- 
termined by standard numerical methods. 

APPENDIX D: ANALYTIC FIRST/LAST CROSSING 
SOLUTIONS IN QUANTIZED, INTERMITTENT 
TURBULENCE 



where Pio(S[S] \B'[S']) is the probability of having S at 5 given 
an earlier value B'[S']; infinite divisibility gives Pio(5[S] | B[S']) = 
Po(S — B'\S — S'). Taking the derivative of both sides with respect 



Here we derive the exact analytic solutions for the first & last- 
crossing distributions in intermittent turbulence as approximated 
in § |7.1| in the text. Consider intermittency of the form character- 
ized by the |She & Waymire| p995) ; |Dubrulle] < [T994] > model, i.e. 
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some characteristic j3 p and 7', with Gr(x) — > S(x) (i.e. pure quan- 
tized log-Poisson processes). Our derivations will closely follow 
those for the first and last-crossing distributions for the Gaussian 
(non-intermittent) cases given in Zhang & Hui ( 2006 ) and Paper II, 
respectively. We refer to those papers for more details. 

Dl First Crossing 

Dl.l Exact Discrete Solution 

Consider first the case with /3 p < 1. Rather than 7', it is more con- 
venient to work in the directly related variable 5 (the variance as 
a function of scale) or A (the variance of the related Poisson dis- 
tribution). We begin at some sufficiently large scale where 5 = 0, 
and consider steps in scale (successively smaller radial averaging) 
that correspond to some increasing 5 — > S + AS. At each scale, we 
can evaluate the density distribution and determine the fraction // 
of trajectories which are experiencing a first crossing in that step. 
Obviously when 5 = 0, 8 = (i.e. all densities are at the mean on 



this scale). From Eq. 



39 



Pi'i = P'p l/V" 2 ) ' Pn> an d fr° m the def- 
inition of A (Eq.|41fr, we have at each scale - equivalently for any 
trajectory - that 



8(R) = In 



(p(R)\ 
V po / 



m(R) In P P + X(R) 
m(R) ln/3 p + S(R) 



-Pp) 

-Pp 



(ln/3 P 



(Dl) 
(D2) 
(D3) 



where m is a Poisson random variable with mean and variance 
A = X(R) = 5/(ln^ p ) 2 (S = 5(6) defined as throughout this pa- 
per as the variance in the logarithmic density field at each scale). 
Because ln/3 p < 0, this value being above some critical density <5 cl i t 
(i.e. crossing the barrier 6(5) requires 



m < m c (X) • 



1 



5(5) -5 



1-0p 
(ln/3 P ) : 



(D4) 



Note here B = In (pan/ po) strictly (there is no ±5/2 term, because 
of how we define the PDF in what follows). 

It must be true that at the "initial" (5 = 0) scale m c < (true 
for any 6 > 0), so that not all trajectories are crossed. A critical 
difference from the Gaussian case is that while 8 can increase or 
decrease in a step, m is a Poisson variable and so is a positive in- 
teger which cannot decrease. So there are no crossings until some 
scale A where m c — > 0; at this point the fraction of all the volume 
with m — (the fraction = exp(— X[m c = 0])) instantly "crosses." 
All trajectories with m > have 8 < B at this scale. Going further 
in scale, A continues to increase. If m c (X) at that scale decreases, it 
is irrelevant, because there can be no decrease in m values for any 
"trajectory," hence anything which now falls above m> m c already 
had a first crossing. If m c (X) eventually increases, it is also irrel- 
evant as long as < m c < 1, because all trajectories which have 
not already crossed (i.e. had m — 0) still lie above m c . However 
at m c = 1 there will suddenly be additional crossings, and so on. 

Crossings can only occur at discrete X mc where m c = 0, 1, 2, 

This means ff is nowhere differentiable. We must define it strictly 
in integral/sum formulation. 

Define the discrete fraction of crossings at each subsequently 
increasing positive integer value of m c , 



AF, 



A [m c 



d5//(5[A]) 



(D5) 



'A[m c -1] 1 

Since crossings occur at the discrete m c , the fraction crossing 



at a given m c must be equal to the expected fraction having m = m c 
at that A, 



P(m I A) = P (m I A) = — exp (-A) 



(D6) 



after subtracting the fraction which previously crossed at m' c < m c 
and now have the value m. For Poisson statistics, the probability to 
have a value m at A, given an earlier value m at A' < A is just 



P(m[X]\m'[X' < A]) = P (m — m | A — A') 
So we obtain the integral equation 



(D7) 



m c -l 



AF f mc = p o(m c I X[m c ]) - AF fkH"ic - 



k=0 



k\X[m c ]-X[k]) 
(D8) 



defined for m c = 1, 2, 3, with 

AF f0 = P {0 1 X[m c = 0]) = exp (-X[m c = 0]) (D9) 

This is defined at the A = X[m c ] with integer values of m c only; 
ff(S) = everywhere else. 

Note that if m c (X) equals a specific integer value at multiple 
5 or A, X[m c ] is defined as the smallest such value of A (i.e. the 
largest physical scale, since we are "counting" in the direction of 
increasing A). 

Once again we note that this function is everywhere non- 
differentiable, therefore //(5) is nowhere defined. However we 
can define the discretely averaged //(5) by assigning a "width" 
A5 corresponding to the AA between m c values, then defining 
(ff(S)) = AF/ m /AS at each m c . For any 6(5), it is now straight- 
forward to evaluate (ff(S)). 

D1.2 Continuum Limit & Linear Barrier Solution 

To compare to our derivation for the Gaussian statistics case, we 
can take the continuum limit of Eq. |D8[ i.e. m> 1 so that the dis- 
creteness effects between m can be ignored. First note that since 
6o(0 1 A — !> 0) — > 1, we can write Eq. 



D8 



^ AF fk P (m c -k\ X[m c ] - X[k]) = P (m c | X[m c ]) (D10) 

<r=0 

From the definition of AFf it is trivial to see that this becomes, in 
the continuum case 

rs 

dS'f f (S')P c (m c [S]-m c [S'] | A — A') = Pc(m c [S] \ A) (Dll) 

where again A = 5/(ln/3 p ) 2 and we now allow m c to assume the 
continuum of values for all 5, and define 



P c (m I A) 



A" 



1 



T[m+l] 



exp 



2A 



(D12) 



where the latter is the limit of the Poisson distribution for large-A. 
Note that in this limit 



6 C K|A) 



/' ; 



(D13) 
(D14) 



|ln/3 p | / (B-jlS) 2 ^ 

■ ' exp — — — 

V2/^S V 25 

l-ff p + ln/3 p 
(ln/3 P ) 2 

i.e. this is just the log-normal probability of 6, for a mean (8) = p, 5, 
and for f3 p — > 1 (no intermittency), ft—i—l /2, exactly as in the case 
of pure lognormal statistics. 

For either form of P c , this is a Volterra integral equation of the 
first kind, and is straightforward to solve via standard numerical 
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methods for any arbitrary barrier 6(5). This is explained in Zhang 
|& Hui|p006) in detail. 

Now consider a linear barrier of the form 6 = In (pcrit/po) = 
Bo + fio S. It is straightforward to verify that the continuum limit 
equation is solved by 

ff(B = \n( Pcdt /po)=B + fioS) = tjP c (m c \X) (D15) 

As P P — > 1, this becomes identical to the exact solution for the log- 
normal (non-intermittent) case, as it should. 

D2 Last Crossing 

D2. 1 Exact Integral-Discrete Solution 

Now consider last-crossings. The derivation is very similar, but now 
we "begin" our evaluation at R — > 0, where 5 takes some finite value 
and B is sufficiently large that all trajectories are below the barrier. 
This means again that m c < must be true (since B>Son these 
scales). Going in the direction of decreasing 5, now, all trajectories 
are un-crossed until m c = 0, when the fraction P = Po (0 1 \[m c = 0] ) 
cross from being "below' the barrier (in terms of 8) on all smaller 
scales to "above" the barrier on a larger scale (i.e. have a last- 
crossing). Note that this is not necessarily the same X[m c — 0] as 
was defined for the first-crossing distribution - that was the small- 
est value of A or 5 (largest spatial scale) where m c = 0, and this 
(because we are evaluating in the opposite direction) is the largest 
such value of A or 5 (smallest spatial scale). 

Unlike in the first-crossing case, because we are evaluating in 
the opposite direction, trajectories can only decrease in m in the di- 
rection we evaluate. But we still are evaluating when trajectories go 
from m > m c Xom< m c . This, now, can occur between integer val- 
ues of m c (e.g. a trajectory could suddenly decrease in m, crossing 
from above to below a non-integer value of m c . 

At some A, it must be true that the sum of all trajectories which 
have "last crossed" from below to above the barrier, plus the sum 
of all trajectories which are currently below m c but have not at any 
point crossed, is unity. So 



1 = 



' fi{S')dS' + n(m|A) (D16) 



where II represents the probability of being at some in > m c , i.e. 
below the barrier, without having previously crossed, so is the stan- 
dard probability Po after subtracting the probability of crossing at a 
previous scale A' and then arriving at the value m 



II(w| A) = Po(m | A) — 



dS'fe{S')Pm (m[A] | FIX( m ;.[A'])) 
(D17) 



Here FIX(m c ) is the integer floor of m c (since an integer crossing 
from above to below m c "lands" at this value). The probability Pqi 
here represents the probability of a transition from m! to m in a step 
in the direction of decreasing A' — > A. This is related to the proba- 
bility Pio of a transition from m to m! in the direction of increasing 
A — > A' by Bayes's theorem 

P m (m[\]\m'[\' > A]) = Pio(m'[\'] | m[X\) ~wy~T\TT\ (D18) 

Po(m' | A') 

Recall Pio(m'[A' > A] j ra[A]) = P (m' - m A' - A) for Poisson 
statistics. Therefore we obtain Poi directly. 

The previous equations can be considerably simplified if we 
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use the following relations: 

FIX(m c ) 

^Po(m|A) = l- ]T Po(m|A) (D19) 

/»>/»£ 111— 

FIX(m c ) 

Poi (m I FIX[ m ;.]) = 1 - Pm ( m I FIX K]) ( D2 °) 



I dS' f e (S')P m (m\FlX[m' c ]) (D21) 
dS'fe(S') Y, fl>i(m|FrX[in£]) 



Expanding II (m \ A) and using the relations above, we obtain 
the following integral equation for ft : 



„\ R=0 FIX(m c ) riA W 

/ dS'MS') J2 Mm[A]|ErX[w£[A']])= p o(«|A) 

« A ,„ — n in — n 



FIX(m r ) 



with 



Poi(m[A]|m'[A']) = 



m'\ A'" (A' -A)" 



m\{m —m\\ 



x ,m> 



(D22) 



(D23) 



Unlike in the first-crossing case, this is differentiable at most 
5, since the sum terms vary continuously in A for fixed FTX(m c ). 
However, when m c crosses a new maximum integer value, FIX(m £ .) 
changes discretely by unity, so both the sum terms change dis- 
cretely. At these values, ft is non-differentiable, and must be de- 
fined only in terms of the integral, i.e. it undergoes a discrete jump 
in the integral AFe as in the first-crossing case. Thus we must again 
only define fe completely over some discrete interval average {ft). 

D2.2 Continuum Limit & Linear Barrier Solution 

In the continuum limit, the sum terms in Eq. |D22| become 

FK(m„) 

2>w^^^[' + »(w)] ,D24) 

m— 



where again the latter is for the A> 1 normal limit of the Poisson 
distribution, and 



FTX(m c ) 



^2 Poi (m[A] | FIXKfA']]) -> - [l +Erf 



m c — (A/A') m'. 
V2(A/A')(A'-A) 
(D25) 



Again in this limit we obtain a differentiable Volterra integral 
equation of the first kind, which can be solved by standard numeri- 
cal methods for any 6(5). 

For a linear barrier again (of the same form as we considered 
for first-crossing), it is straightforward, albeit tedious, to verify (by 
taking the derivative of both sides of this equation and solving the 
resulting Volterra integral equation of the second kind) that the con- 
tinuum solution is given by 

/f(6 = ln(p qrit /po) =6 + ^o5) ^ \fXo+p,\P c (m c \X) (D26) 

Here note that since the true "run" in the numerator of v is given 
by v = (6o + (/io + p) 5) /5 1 ' 2 , this is the same as the form of the 
solution for lognormal (non-intermittent) statistics, but with an ap- 
propriate modification jl for the /3 p -dependent run in the mean. In 
the limit P P —¥\, this becomes identical to the solution for the non- 
intermittent case. 
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D3 Inverse (3 Models 

For P P > 1, the relation between S and m (in terms of f3 p and A) is 
identical, and we can define the same m c . The difference is, since 
ln/3 p > 0, 5 now increases with m. So 5 > B means m > m c . This 
qualitatively changes our previous solutions. 

For first-crossing, we again begin at m = and 5 = 0; now 
B > is required so the mass is not already crossed, but this sim- 
ply means m c > 0. As we go down in scale, m can again only 
monotonically increase; but since m > m c gives 8 > B, crossings 
can occur "between" integer values of m c . This is symmetric to the 
last-crossing distribution with R p < 1, though still "stepping" in the 
direction of increasing A (using Pio instead of Pqi), and switching 
the sum over II from m > m c to m < m c . This gives: 

/ dS'f e (S') V P (m-FIX + [m c ] A- A') = V P (m\ A) 

Jo >, ^ 

m>m c m>m c 

(D27) 



m>m c 



where we use FiX + (m£) to denote the "rounded up" value of m c , 
i.e. the minimum integer value n > such that n > m c . 

For last-crossing, this reversal makes the situation symmetric 
to the f3 p < 1 first-crossing solution, and is nowhere differentiable. 
This gives 

AFt mc =P [m c \\)- AF ek P Q1 {m c [\]\k[X k ]) (D28) 

/,=/«, +i 



APPENDIX E: THE LOG-DENSITY VARIANCE AS A 
FUNCTION OF SCALE 



In Eq. 12 we approximate the run of density variance S(R) with 
scale by 



S(R, p) 



\W(k,R)\ 2 In [ 



b 2 vj(k) 



s (p,k) 2 + >t 2 k- 



dlnfc (El) 



Testing this approximation in numerical simulations (specifically, 
how to map between turbulent velocity power spectrum and den- 
sity variance as a function of scale) is an important question for 
future study. Subsequent work should study this approximation not 
just in simple isothermal cases, but also its extensions to vertically- 
stratified, rotating shear flows with non-isothermal gas. 

That said, we can make some very preliminary comparisons 
with published simulations. Most such studies focus on idealized 
driven turbulent boxes, with isothermal gas and no rotation. In this 
case, Eq,|12|becomes 



S(R) = 



\W(k,R)\ 2 In 



(E2) 



where k max corresponds to the box scale. Without loss of gen- 
erality, we can use units where k milx — 1. Note, by definition, 
bvt/cs — M compressive, i.e. only the compressive component of the 
velocity enters. So we can simplify by stating all quantities in terms 

of Mi compressive I 



and Ai with -Mcompressive ~ b Ai, as discussed in §|5] Determining the most 
accurate mapping between the two is another important question for future 
study, but for clarity here, we separate it from the distinct question of the 
mapping between M compressive 

and S(R). 



By definition, S(R) is just the convolution over the log- 
density power spectrum with the window function (S(R) 
J\W(k)\ 2 S lnp (k)dk). SoEq. 



E2 



5*111 p (k) — k 
= k~ 



ln[l 
ln[l 



implies a power spectrum 

M compressive \ 

M;, k , - p ] 



(E3) 



where the latter assumes a power-law velocity spectrum and 

M\.c. — ^^compressive (k — ^max). 

Eq. |E3| follows from the assumption that a Mach number- 
dispersion relation 5 



\n[l+ Mi 



Konstandin et al. 



Compressive] (^ - § 

|2012b| >, applies scale-by-scale. This follows if we assume the den 
sity field is strictly self-similar (see Paper I for details). Although 
reasonable, this is not unique. 

For example, consider instead a model where S(R — > 0) ~ 
ln[l + Ml o] is true exactly, but only as a box-integrated state- 
ment about the variance on small scales. Using S(R — » 0) = 
S]„ p (k)dk, and M 2 C ^ = f™E v , c (k)/c 2 dk (here E v . c is the 
power spectrum of compressive velocity fluctuations), we estimate 



fait /,i (p- l)kE vc c~ s 

Jlnp W 



Ml k- 



k(l+kE v , c ci 2 ) l+Miok'-i'ip-l)- 



(E4) 



In the sub-sonic regime (M c ,o < 1 or ^> 1), all of these 
scalings become identical. In either case, §i up (k) — > M 2 cA k~ p oc 
E VlC (k), the velocity power spectrum. Since the amplitude of fluctu- 
ations is small, the density and log-density power spectra are equiv- 
alent: S p (k) oc S\n P (k) oc E VtC (k) oc k~ p . This is a generic and well- 
known prediction in the weakly-compressible regime (see |MoriF] 
|gomery et al.|1987 1, which has been verified in a range of simula- 
tions (see e.g. Kowal et al.[ 2007 , Schmidt et al. 20091. The same 
conclusion should generalize as predicted here for polytropic equa- 
tions of state in the weakly compressible limit (Biskamp 20031. 
And the box-integrated variance S(R — > 0) — > ln[l +Mi ] ~ 
Ml o, in agreement with the standard Mach number- variance re- 
lation. 

The most appropriate scaling when M(k) S> 1 is more am- 
biguous. Federrafh et al. (2010) directly measure E v , c (k) and 
Si Bp (k) in simulations of compressively and solenoidally-driven 
turbulence with M ~ 5.5. We can therefore directly compare 
to the measured S\„ p (k), using the m easu r ed E v , c (k) 
e (k) = c7 2 j k °° E V)C (k) dk). Both Eq. 



Eqs. 



E3IE4 



E3 



& 



E4 



agree 



(and compressive V 

reasonably well with both compressive and solenoidal simulations 
(within ~ 20% at k outside a factor « 2 of the driving and reso- 
lution scale). Interestingly, Eq. |E3| agrees slightly better with the 
compressive case, while Eq. |E4| agrees slightly better with the 
solenoidal case. We can perform a similar exercise with the sim- 
ulations in |Kowal et 10^2007) , spanning M ~ 0.2 — 7 and Alfven 
Mach number Ma ~ 0.7 — 7; here, however, E r , c (k) is not mea- 
sured (only E v (k), the total velocity power spe ctrum), so we ap- 
proximate E v , c (k) ~ b 2 E v (k). In this case, Eqs. E3|E4 are consis- 
tent with the simulated S\ np within a factor ~ 2 (so much of this 
deviation may owe to b not being constant). 

Note that, by construction, Eq. |E4| reproduces the "usual" 
box-integrated variance-Mach number relation S(R — > 0) 



E3 



ln[l +M C oL in both sub and super-sonic limits. However Eq. 
leads to a slightly different box-integrated relation when M c ,o 2> 1 
for a power-law E v , c (k), this is the dilogarithm. This is identical at 
sub-sonic M c .o but rises more steeply (albeit still logarithmically) 
at large M c ,o- In Hopkins (2012a), we find (from a large compi- 
lation of simulations) that the mass-weighted statistics appear to 
more accurately follow the logarithmic scaling, while the volume- 
weighted statistics (which we use in this paper) more closely follow 
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the dilogarithm scaling. This motivates our particular choice in this 
paper. However, the simulation results clearly remain ambiguous, 
and their disagreement may be related to more fundamental aspects 
of the density PDF, such as intermittency. 

In either case, these scalings are broadly consistent with obser- 
vations of the projected surface density power spectrum in galactic 
gas JStanimirovic et al.|1999||Padoan et al.|2TjQ6||Block et al.|2010| >. 
However, with the present uncertainties, the ability to distinguish 
between the two is limited. 

Given this uncertainty, we have re-run our key calculations 
in this paper with Eq. |E4| (re-adding the appropriate k terms), in- 
stead of our usual Eq. |E2| We find no change in any of our qualita- 
tive conclusions. Quantitatively, the correlation functions and last- 
crossing distributions are largely unchanged; first-crossing distri- 
butions tend to "shift" slightly in mass, but at the factor < 2 level. 

APPENDIX F: ON THE CONVOLUTION OF 
LOGNORMALS 

In this paper, we have approximated the real-space density distri- 
bution as a lognormal averaged on all radial scales. But this can- 
not be exact if the real-space fluctuations on each scale are un- 
correlated. In real-space, mass conservation means that the aver- 
age density in the summed volume of two sub-regions must be 
(V\ +Vi)pi+2 = V[ pi + Vi pi (with Vi , V% the volume of respec- 
tive sub-regions); i.e. the summed real-space variable is the linear 
p, not the logarithm In p. 

But it is well-known from numerical calculations that the log- 
normal approximation is quite accurate (especially for the high- 
density tail of the density distribution, our focus in this paper), 
over a wide range of spatial and resolution scales (including re- 
averaging the same simulations over smoothing scales; see e.g. 



Vazquez-Semadeni||1994||Padoan et al. |1997 


Scalo et al. 


1998 


Nordlund & Padoan|1999||Ostriker et al.|1999| 


Kowal et al. 


|2007 



calculating the density PDF at different resolution would yield fun- 
damentally different shapes. 

The key is the well-known mathematical result that the con- 
volution of two lognormal distributions (i.e. the sum of linear vari- 
ables drawn from lognormals) can be very well-approximated by 
a lognormal (related to their similar characteristic functions). In 
particular, the Fenton- Wilkinson approximation lFent on|1960} , in 
which the convolution of lognormals is approximated as a lognor- 
mal, captures the high-density tail (positive fluctuations) very ac- 
curately (to within w 2% at p > (p) ; see e.g. |Mehta et al.|2007) 
given the simple constraint that the first two linear moments of the 
approximate lognormal are matched to the exact convolution result. 
If we satisfy these moment conditions, then, the residual deviations 
from log-normal (owing to the convolution) are much smaller than 
those owing to realistic levels of intermittency or small deviations 
from 7 = 1, let alone dynamic effects of collapsing regions, self- 
gravity, and large-scale perturbations (global modes) not captured 
in our analysis. 

For the first moment, our enforcement of mass conservation 
in the "steps" of the lognormal (by setting the median value in the 
PDF to —5/2 for a lognormal) means that this always automatically 
satisfied in real-space. The second moment condition implies (in 
Fourier space) 

<(ln(p/po) - <ln (p/po)» 2 ) = ln(l + {(p/p - <p/po» 2 » 

i.e. S[lnp] = ln(l + S[p]) (where S[lnp], S[p] are the variance in 
Inp and the linear p, respectively). But this is used to derive 5[lnp] 



itself (see e.g. Paper I, Padoan & Nordlund 2002, Konstandin et al. 
|2012b| and references therein), so is also automatically satisfied 
to the accuracy that our Eq. [12] represents the log-density power 
spectrum. 

Moreover, the non-lognormal intermittency models in the text 
exactly obey the correct linear convolution condition, for specific 



choices of 3 and 7' or T in Appendix [B| (see Appendix A in Hop- 
|kins|2012a> . 

If, however, the behavior of the low-density tail of the density 
distribution is desired (for example, calculating the mass function 
of underdense "voids" or ISM "holes," as in Paper I), this approxi- 
mation becomes much less accurate. This is confirmed in turbulent 
box simulations (Federrath et al. 2010 1. In this limit, one should in- 
stead adopt the Schwartz- Yeh (Schwartz & Yeh 19821 approxima- 
tion (matching the first two moments in log-space), which becomes 
similarly accurate for p <C (p). 

If the volumes summed are discrete real-space spheres (as op- 
posed to our default choice of a the Fourier-space top-hat), some 
additional care is needed, discussed below (Appendix [G|. 

APPENDIX G: THE WINDOW FUNCTION & 
REAL-SPACE AVERAGING 

In § [2] and throughout this paper, we simplify by assuming that 
densities are averaged over a window function which is a top-hat in 
Fourier (k) space (W(k, I ty in Eq .|l2| = 1 for k<R-' and = for 



Bond et al. 



1991b, this is what allows us 



k > R~ ). As discussed in 
to treat fluctuations between scales as an uncorrelated random walk 
(at least in the isothermal case). This is necessary to obtain closed- 
form analytic solutions for many of the quantities we consider. 

However, in certain situations it might be more appropriate 
to use a different filter function, for example a simple real-space 
sphere (a real-space top-hat of radius R). But this necessarily means 
fluctuations between scales are correlated, since "incrementing" the 
averaging scale integrates over contributions from all fc-modes|^] 
Mathematically, this is related to the Fourier transform of the win- 
dow function. The overdensity in real-space is defined by: 

(lnp(x,^,„)) = / dx'W(\x' -x\,R w ) lnp(x') (Gl) 



so in fc-space lnp(k, R w ) = W(k, R w ) lnp(k). For the spherical vol- 
ume in real-space, W(x = |x' — x|, R w ) — (47r^,./3) _l for* < R w 
and = for x > R„, which has Fourier transform 



W(k=\k\,R H .) 



3 [sin (kR w ) - kR w cos (kR w )] 
(kR~J 3 



(G2) 



31 Consider the sum of a specific pair of discrete real-space volumes 
V\ and V%, which have approximately lognormal PDFs for each indi- 
vidual density within the volume, into a larger volume V1+2 = V*i + V*2 
with mean density P1+2, and consider the PDF of density that results as 
an approximate lognormal (with the Fenton- Wilkinson approximation in 
§ [Fj. The first moment condition is just a restatement of mass conserva- 
tion: V"i+2(pi+2) = V\ (p\) + V*2(P2}- By our definition of the median 
—5/2 of each lognormal "step", (pi) = (P2} = (pt+2)> so this is al- 
ways satisfied. Matching the second moments in the mass conservation- 
restricted sum, and simplifying, we obtain S[lnpi+2] = ln(l+Sj where 
§ = fl (exp[5i] - 1) + /| (exp[S 2 ] - 1) - 2/i/ 2 ({piPl/pD - 1) with 
fi = Vj/Vi+2 and 5, = 5[lnp,]. It is then straightforward to show that, if 
the variables pi and P2 are uncorrelated, there is only one possible power 
spectrum shape for S, namely that of integrated Poisson fluctuations (with 
fluctuation "number" proportional to volume). The existence of any non- 
trivial structure in the power spectrum therefore necessarily implies corre- 
lated steps in real-space. 
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Figure Fl. Mass functions, as Fig. [3] with three different window func- 
tions W(x, R) used to compute the real-space density fluctuations on a scale 
R (see Appendix |Gj). We compare a sharp Fourier (&)-space top-hat (the 
default choice in this paper, because it allows exact analytic solutions for 
the MF), a Gaussian filter (Gaussian weighting in real-space R and Fourier- 
space k), and a real-space top-hat (sharp-edged sphere in real-space). Top: 
MF for our standard model with these window functions, with no correc- 
tions applied to any quantity. The MFs are nearly identical in shape, but off- 
set in normalization/mass scale at high-mass - this is because the definition 
of mass & volume for an "object" are different. Bottom: Same, but enforc- 
ing the same definition of mass (& mass- variance relation) for each window 
function. The larger normalization offsets go away. The Gaussian and real- 
space sphere windows (themselves giving nearly identical results) produce 
smoother trajectories in Fig. [TJ (feature less small-scale noise/structure), so 
produce fewer "back-and-forth" crossings at intermediate scales. This sup- 
presses the low-mass end of the first-crossing MF and high-mass end of the 
last-crossing MF. 



Thus calculating the increment in In p in real-space involves inte- 
grating over contributions from all modes in £-space, and "steps" 
are correlated. 

The introduction of a non-trivial window function changes our 
calculation in two ways. First, we must modify the calculation of 
S(R) in Eq. 



12 



inserting the appropriate W (k, R). Second, we must 
modify our calculation of the density field. Derivations of the ap- 
propriate walk are given in |Bond et al.| l fT99T) ;[Zentner (2007]l, but 
the result is simple. For the Monte-Carlo method of evaluating "tra- 
jectories" in the (zero-mean) variable 5 (R), each trajectory obeys 
the integrated Langevin equation: 



s(r) = r 

Jo 



TZ(\nk)W(k,R)dlnk 



(G3) 



where 1Z is the "stochastic force," a Gaussian random variable 
(independently drawn at each dlnfe interval) with zero mean and 
variance (TL 2 {\nk)) = (dS / din k)/\W(k, R)\ 2 (i.e. the integrand in 
Eq.[T2] without the window function). For a Monte carlo ensemble, 



each trajectory S(R) should be calculated with its own independent 
set of 7?.(lnA:) for all R (but that set is preserved for each R used in 
the integration) ^] 

In Fig. [FT we re-calculate the mass functions for our "stan- 
dard" model, using three different window functions: the £-space 
top-hat (standard in the text), the sphere (real-space top-hat) above, 
and an intermediate (Gaussian) window function]^] First, we sim- 
ply insert these forms of W into the appropriate equations for 5 and 
S, but keep the rest of the model fixed. As shown, this leads to pre- 
dicted MFs with very similar shape. Interestingly, the Gaussian and 
real-space sphere filters predict nearly identical MFs to one another, 
but with some differences from the sharp &-space top-hat. The last- 
crossing MF predicted in the former cases has a slightly steeper 
slope (but recall the very large dynamic range plotted and note the 
difference is quite small, from a logarithmic slope of ~ —2.0 to 
~ —2.1). And the first-crossing MF, while nearly identical in shape, 
is offset in the former cases to higher masses (by an apparently 
substantial factor of ~ 3). However, as cautioned by Bon d et al.| 
{199TJ and[Zentner (2007), care is needed in this simple compari- 
son, because for otherwise fixed properties, the definition of "effec- 
tive mass" and volume associated with each of these filters is in fact 
different, as is (as a consequence) the relation between mass or vol- 
ume scale and the variance associated with that scale. Therefore it 
may be more appropriate to compare the different filters, enforcing 
a choice of filter radius such that the mass-volume, and volume- 
variance (or mass-variance) relations are the same. If we do this, 
also shown in Fig. |Fl| the large normalization difference, and some 
of the difference in slope, disappears. 

There is still a non-trivial remaining difference, in the sense 
that the last-crossing MF has a slightly steeper slope with a Gaus- 
sian or real-space sphere filter (less mass at the highest masses) and 
the first-crossing MF has a slope modified in the opposite sense 
(less mass at the lowest masses), but this is the regime which does 
not contain most of the mass in either MF. What happens here is 
that these filters, by smoothing the contribution of Fourier modes 
from all scales, exhibit less small-scale "structure" - the trajecto- 
ries 5(R) are considerably more smooth than those in Fig.[TJ while 
exhibiting identical variance. As a result, there is less probability 
of intermediate scale, large random fluctuations which are con- 
centrated over a narrow range in spatial scale R (i.e. have a first- 
crossing just slightly larger in scale than last-crossing). In other 
words, these different filters further enhance the tendency discussed 
in the text for the dynamic range of fragmentation to be concen- 
trated in first-crossings at the maximal instability scale and in last- 
crossings at the sonic scale. 

One potentially important caveat to this comparison is that we 
have not modified the collapse threshold/barrier (i.e. Eq,|14} for the 
different window functions. Ideally, one would re-derive this for the 
appropriate window function. The form we use is based on the dis- 
persion relation derived for single Fourier (k) modes, so it is most 
appropriate to link this to a fc-space top-hat window function (as 
we have done in the text). It would be interesting to see whether re- 
deriving this for Gaussian or real-space sphere overdensities would 
increase or minimize some of the differences seen here. Moreover, 



32 For a more detailed discussion of how to treat even more complicated 
window functions via the path-integral formulation, see Maggiore & Riotto 
gOTOa). 



This is often used in calculating 
and is also numerically convenient. It 
exp(-x 2 /24)/([27r]V24) and W(k, R w ) 



"smoothed" density fields, 
is defined by W(x, R w ) = 



--enp{-k 2 /2R~ 2 ). 
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although windows such as a real-space sphere seem physically sen- 
sible, in an ideal case the window function (and barrier derivation) 
would be matched to some typical "shape" of density fluctuations, 
which is certainly not obvious and may not be universal (since we 
are considering the fully non-linear density field). 



Gravo -Turbulent Fragmentation 
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